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I .  INTRODUCTION 


This  report  presents  and  discusses  a  general-purpose  FORTRAN 
equation-fitting  program  called  FINLIE. 

Assume  that  the  behavior  of  some  physical  system  can  be  adequately 
described  by  a  set  of  equations  involving  one  independent  variable  x 
and  N2  dependent  variables  (N2  >  1).  FINLIE  requires  that  these  equa¬ 
tions  be  reducible  to  one  of  two  forms: 

Ca)  a  system  of  N2  first-order  ordinary  differential  equations 
of  the  form 

dy./dx  =  f.  (X,  Y.  C)  [j=l.  2,  ...  N2]  (1) 

where  Y  is  the  vector  of  N2  dependent  variaL  es: 

Y  =  y2»  •••  yjs}2^ 

and  where  C  is  a  vector  of  N3  linearly  independent  parameters  (N3  >  0) : 

C  H  (Cj,  c^,  ... 

(b)  a  system  of  N2  algebraic  and/or  transcendental  equations  of 
the  form 


yj  =  gj  (X.  Yq,  C)  [j  =  1,  2,  ...  N2]  C2) 

where  Y^  is  the  initial  condition  vector: 

720'  •  •  •  Yi^2,0^ 

The  user  writes  his  system  (1)  or  (2)  as  a  FORTRAN  subroutine  whose 
name  is  submitted  to  FINLIE.  FINLIE’ s  task  is  to  adjust  the  parameters 
and  initial  conditions  of  (1)  or  (2)  so  as  to  fit  the  solution  curves 
to  measurements  taken  on  one  or  more  of  the  dependent  variables.  For 
system  (1),  no  knowledge  of  the  form  of  the  solution  is  necessary. 
Indeed,  we  may  in  general  assume  that  system  (1)  possesses  no 
closed-form  solution  of  the  form  (2).  Otherwise,  we  would  fit  the 
solution  equations  rather  than  the  differential  equations. 

System  (1)  can  be  linear  or  nonlinear  in  the  parameters;  system  (2) 
can  be  linear  or  nonlinear  in  the  parameters  and  in  the  initial  condi¬ 
tions.  However,  linear  parameters  and  initial  conditions  are  not  much 
of  a  challenge  to  FINLIE.  Indeed,  the  word  FINLIE  can  be  viewed  as  an 
acronym  for  "Fitting  NonLInear  Equations";  the  program  was  created  to 
handle  nonlinear  situations.  (System  (1)  may  also  be  nonlinear  in  the 
more  common  sense  of  "nonlinear  in  the  dependent  variables";  for  our 
purposes,  this  is  irrelevant.) 
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As  a  rather  elementary  example  of  system  (1),  consider: 


dy^/dx  =  l/y2 

dy^/dx  =  -Cj  ^  ^ 


C3) 


Here  N2  =  N3  »  2.  If  x  is  interpreted  as  distance,  y^  as  time  and 

y^  as  the  magnitude  of  a  missile's  velocity,  then  (3)  is  essentially 

the  drag  equation  for  a  horizontal  flight  in  which  the  drag  coefficient 
varies  linearly  with  Mach  number. 


One  of  the  reasons  we  chose  this  particular  example  is  that  it 
does  possess  a  closed-form  solution: 


^1  "  ^10  ■  ^^2 
y^  =  (bu  -  C2)"^ 

where  u  =  exp  [c^  (x  -  x^)] 

b  =  =2  ^  (720^'^ 


(4) 


In  "real  life"  we  would  always  fit  (4) --which  is  of  the  form  (2)-- 
and  forget  about  (3).  In  this  report,  however,  we  will  use  both  C3) 
and  (4)  to  illustrate  our  remarks. 


FINLIE  is  given  measurements  on  the  first  N1  of  the  N2  dependent 
variables;  that  is,  on  y^,  y^  ...  y^^j  where  1  <  N1  <  N2.  The  m-th 

data  point  thus  consists  of  N1  measurements  at  the  independent 

variable  value  x  : 

m 

^Im’  ^2m’  ^Nl,m^ 

where  y.  denotes  the  measured  value  of  y.  at  x  . 

■'jm  •'j  m 

Assume  that  the  measurements  have  been  obtained  from  one  or  more 
distinct  experiments,  each  experiment  having  its  own  initial  condition 
vector.  Because  our  first  practical  application  of  FINLIE  was  to  rounds 
fired  in  an  enclosed  range,  we  will  call  each  distinct  experiment  a 
round.  By  "multi-round"  data,  then,  we  mean  NR  sets  of  measurements 
(nr  ^  1),  all  applicable  to  the  same  system  of  equations  and  hence  help¬ 
ing  to  determine  the  single  parameter  vector  C,  but  each  measurement  set 
determining  its  own  initial  condition  vector. 


Thus  there  are  NR  x  N2  initial  conditions  to  be  determined: 


IC  =  {(Y^)j,  (Y^)2. 
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FINLIE  requires  that  these  initial  conditions  refer  to  the  same  independ¬ 
ent  variable  value  for  every  rcund.  However,  need  not  coincide 

with  any  value  x  at  which  measurements  were  taken  and  x,,  need  not  even 
'  m  o 

fall  within  the  interval  bounded  by  the  smallest  and  largest  of  the 
values.  (Of  course,  the  farther  lies  from  that  inteirval,  the  more 
unreliable  is  the  extrapolation  to  that  point.) 

Ws  assume  that  within  each  round,  the  x  values  increase  with 

m 

increasing  m.  For  example,  if  we  have  two  rounds  with  four  and  five 
data  points,  respectively,  then 

X,  <  x-  <  x,  <  x, 

12  3  4 

and 

5  6  7  8  9 

but  no  demanas  are  made  on  the  combined  ordering  of  the  nine  values. 

A  member  of  the  first  string  of  inequalities  above  can  be  less  than, 
equal  to  or  greater  than  some  member  of  the  second  string. 

For  convenience  we  coin  the  word  "paramic"  to  mean  "parameter  or 
initial  condition."  Of  course,  the  initial  conditions  are  parameters 
of  a  sort:  parameters  whose  values  can  change  with  x^  and  with  the 

round.  Thus,  for  example,  system  (4)  could  have  been  written  in  terms 
of  four  "parameters";  say,  in  the  form 

yi  =  Cj  -  CjX  +  (C^/Cj)  z 

where  z  =  exp  f^Cj^x)  ,  This  form  conceals  the  fact  that  the  values  of 
two  of  the  four  c^. '  s  will  change  with  the  initial  conditions. 

By  our  definition,  a  parameter  is  independent  cf  the  choice  of 
Xq  and  applies  to  (and  is  influenced  by  the  measurements  from)  all  the 

rounds.  This  is  the  essential  condition  we  impose  on  the  NR  rounds  to 
be  fitted  simultaneously:  that  the  same  parameter  vector  C  applies  to 
each  round.  The  measured  data  for  any  one  round  may  be  incapable  of 
determining  C  adequately;  the  combined  rounds  have  a  much  better 
chance . 

FINLIE' s  task  is  to  find  the  set  of  paramics 

P  =  (IC,  Cl  (S') 

that  best  fits  the  solution  curver  to  the  multi-round  measurements. 


Note  that  P  consists  of  NR  x N2  initial  conditions  and  N3  parameters,  a 
total  of 


N  =  (NR  X  N2)  +  N3 


(6) 


paramics.  By  a  "best  fit",  we  mean  a  least  squares  fit.  That  is, 
FINLIE  seeks  a  particular  set  P--call  it  ?-~that  minimizes  e,  the  sum 
of  the  weighted  squares  of  the  residuals  of  the  fit: 


N4  N1  _  ~ 

f. (P)  \  ''  w.  [y.  -  y.  (x  ,  P)] 

^  '  jra  ^^jm  ■'j  m’ 


(7) 


where 


N4  =  the  total  number  of  data  points  for  all  the  rounds; 
w^^  =  a  non-negative  weighting  factor  associated  with 


y.  evaluated  at  x^,  using  the  current  value  of  P. 


Other  convenient  measures  of  the  goodness  of  fit  include; 

2  .  cCP) 

(a")  the  estimated  variance  of  the  fit  =  s  = 

(b)  the  estimated  standard  deviation  of  the  fit  =  s 

(c)  the  estimated  probable  error  of  the  fit  =  0.67449  s. 

Note  that  for  a  least  squares  fit  we  must  have  N4  >  N;  that  is, 

there  must  be  more  data  points  than  paramics.  (We  also  assume  that  the 
number  of  data  points  in  each  round  exceeds  N2,  the  number  of  initial 
conditions  for  each  round.) 


The  function  e  is  nondimensional.  Hence,  if  we  let 

[  5  dimensions  of  [  ] , 

Eq .  (7)  implies  that 

fw.  =  (y-'‘L 

'  jm-'d  M 

If  the  user  fails  to  specify  the  values  of  the  weights,  FTNI IE  will  set 
all  weights  to  unity.  This  may  or  may  not  be  adequate.  Usually  the 
weights  are  chosen  so  that  each  term  in  (7)  is  of  the  same  order  of 
magnitude.  This  can  be  done  by  making  w^^  inversely  proportional  to 

the  square  of  the  uncertainty  in  measurement  yjjj^'- 
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where  K  is  a  nondimensional,  positive- -but  otherwise  arbitrary--number . 
That  is,  in  general  only  relative  uncertainties  are  nredea."  suppose, 
for  example,  that  there  are  two  measured  variables: 


y,  (furlongs) ,  for  which  the  uncertainty  in  each  measurement 
^  is  about  ten  furlongs; 

(fortnights),  for  which  each  uncertainty  is  about  0,1  fortnight. 

If  we  choose  K  equal  to,  say,  (a^^)2  in  (9),  we  have 

w,  =  100/100  =  1  (furlong)  ^ 

Im 

w„  =  100/0.01  =  lo'^  (fortnight)"^ 
zm 

Any  other  weights  for  which  '^2m''^'"lm  ~  would  work  as  well.  In  fact, 

any  weights  for  which  the  ratio  is  "close"  to  lo'^--say,  within  a  factor 
of  two  larger  or  smaller — would  probably  work  as  well.  Letting  FINLIE 
set  all  weights  at  unity,  on  the  other  hand,  would  not  work  well  at  all 
in  this  situation.  The  y^  measurements  would  then  have  much  too  great 

an  influence  on  the  fit;  their  noise  would  drown  out  the  y2  measurements. 


If  measurements  are  taken  on  more  than  one  dependent  variable 
(that  is,  if  N1>1),  it  may  happen  that  for  some  data  point  R^,  one  or 

more  (but  not  all)  of  the  measurements  is  missing  or  is  clearly  very 
wrong.  There  is  no  need  to  discard  the  entire  data  point;  it  suffices 
to  set  the  weights  of  any  missing  or  outlier  measurements  at  zero. 


If  we  are  fitting  the  solution  system  (2)  to  the  data,  FINLIF. 
computes  the  values  TjC^jj^.F)  in  (7)  directly  from  the  given  expressions. 

If  we  are  fitting  the  differential  equation  system  (1),  however,  then 
FINLIE  must  obtain  y^Cx^.P)  by  numerical  integration.  When  we  have 

a  choice,  we  pick  (2)  over  (1)  to  avoiiJ  this  integration:  'tis  a 
summation  devoutly  to  be  missed. 


Each  time  FINLIF.  is  called  by  the  user,  it  performs  one  iteration 
of  its  search  procedure.  That  is,  the  user  gives  FINLIE  the  paramic 
set  and  FINLIE  returns  a  set  Pj.  Pj  is  almost  certainly  not  the 

desired  solution,  but  it  should  be  an  improvement  over  Pj^^  in  the  sense 

that  c{Pj)  <  ^(Pq).  The  user  then  gives  FINLIE  the  set  Pj  and  gets 

back  P-,,  and  so  on.  The  process  stops  when  a  specified  convergence 
criterion  is  satisfied  or  some  computational  disaster  arises. 


^HoueVt-v,  for  an  abaa !  uta  fn.tt'i'rnctat  fan 
t>aiseJ  on  c)  K  shauU  ,'<■  1. 


(\ma  aa:  < 


vrof  ntoi'uf,  .' 
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To  illustrate  some  of  the  above  generalities,  we  I'etum  to  our 
sample  systems  (3)  and  (4).  Suppose  th^t  from  three  enclosed-range 
firings  we  obtain  the  data  points  (x^,  listed  in  Table  I. 

Assume  that  the  x  values  in  the  table  are  exact  but  that  each  of  the 
m 

sixteen  y^  measurements  has  an  associated  uncertainty  (seconds) . 


Table  I.  Sample  Data 

Points  for  System  (3)  or 

(4) 

m 

X  (metres) 
m 

yj^(seconds) 

1 

0.0 

2.0000000 

2 

1.0 

2.0100507  I 

3 

2.0 

2.0202034 

Round  El 

4 

3.0 

2.0304591  1 

5 

4.0 

2.0408189 

1 

6 

-3.0 

-0.0147728  1 

7 

-2.0 

-0.0098987  1 

8 

-1.0 

-0.0049746  1 

f  Round  E2 

9 

0.5 

0.0025064 

10 

1.5 

0.007557?  1 

11 

2.0 

0.0101027 

12 

0.  0 

3.0000000 

13 

1.0 

3.0033506 

14 

2.0 

3.0067358 

15 

3.0 

3.0701561 

j  Round  E3 

16 

5.0 

3.0171031  1 

Here  we  have  NR  =  3  rounds 

(the  three  firings),  N4  = 

16  measure 

ments  and  N  =  8  paramics.  The  paramics  are  the  six  initial  conditions 
and  two  parameters: 


^  ^yi0’^20^El 


*'^10’^20^E2’ 


(y 


10’^20^E3’ 


c,} 


(10) 


where  we  arbitrarily  let  x  --the  x  value  at  which  all  six  initial 

o 

conditions  apply--be  zero.  The  values  of  the  eight  paramics  arc  to 
be  adjusted  so  as  to  minimize 

_  1 

w,  [y,  -  y,  (x  ,  P)  ]  " 

Im  Im  1  m  ‘ 


Whenever  only  one  dependent  variable  has  been  measured  (N1  =  1), 
tne  user--unless  he  has  information  to  the  contrary--can  assume  that 
all  the  uncertainties  a'ce  equal.  This  simplifies  matters  by 

allowing  the  user  to  set  w,  =1  for  all  m.  Thus,  for  Table  I,  we  set 
*  Im  ’ 
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_2 

=  1  (seconds)  [m-1, 2, . . . 16] 

The  "measured"  values  in  Table  I  were  actually  obtained  by 

rounding  to  seven  decimal  places  the  values  computed  from  the  solution 
system  (4),  using  Xq  =  0  and 

P  -  {(2,100)gj.  (0,200)p2.  (3, 300)^2*  0.01,  0.0001) 

The  y^^^  values  in  the  table  are  thus  equal  to  number 

of  decimal  places  shown.  FINLIC's  task--given  system  (3)  or  (4)  r’.nd 
the  Table  I  data~-would  be  to  find  P. 

FINLIE  must  be  given  another  bit  of  information  before  it  can 
begin  its  search  for  P;  a  starting  point  P^.  For  systems  (3)  and 

(4)  and  the  Table  I  data,  we  gave  FINLIE  the  relatively  poor  first 
estimate 


Pq  =  {(1.5,50)gj,  (-0.5,250)p2.  f2. 5, 250)^5,  0.02,  0} 

FINLIE  then  proceeded  from  to  P^  to  P2  and  so  on  to  P^,  an  acceptable 

approximation  to  P  (see  Table  II).  Within  the  idiosyncracies  of  machine 
computation,  this  path  from  P^  to  P^  is  the  same  whether  we  fit  system 

(3)  or  system  (4) .  As  one  might  expect  in  a  convergent  situation,  the 
last  two  points  (P^  and  P^)  are  practically  coincident.  Tlie  slight 

discrepancy  between  P^  and  P  is  due  almost  entirely  to  the  round-off 

error  in  the  y^^^  data  of  Table  I. 

Unfortunately,  a  poor  choice  of  Pq  can  sometimes  prevent 

FINLIE' s  ever  finding  P.  Hence  a  reasonable  amount  of  labor  expended 
in  determining  P^  may  pay  dividends.  For  frequently  recurring 

ajiplications,  it  may  be  worthwhile  for  the  user  to  write  his  own 
FORTRAN  subroutine  for  extracting  a  first  estimate  from  the  data 

points.  Usually  only  a  few  of  the  paramic  estimates  are  critical 
for  obtaining  convergence  to  P;  the  remaining  paramics  can  have  sur - 
orisingly  poor  first  estimates  with  impunity.  And  for  some  systems 
uf  equations,  the  choice  of  P  is  very  nearly  immaterial:  all  roads 
lead  to  P. 

A  useful  feature  of  FINLIE  is  its  ability--at  the  user's 
request — to  hold  fixed  the  input  values  of  any  specified  paramics, 
rather  than  allow  those  input  values  to  be  adjusted  by  the  fitting 
]  rocess.  Thus,  for  example,  the  effect  of  a  given  parameter- -say, 
c,  in  system  (3)  or  (4)--can  be  suppressed  during  a  computer  run  by 

giving  that  parameter  an  initial  value  of  zero  and  specifying  that  this 
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Table  II.  Path  from  to  P.^  for  System  (3)  or  (4)  and 
the  Data  of  Table  I. 


yjoCs) 


Round  El _ 

720  (m/s) 


Round  E2 


1.5 

1.998 

1.9999830 

1.9999996 

2.0000001 

2 

2 

2 


50 

69.75 

90.83 

99.18 

99.993 

99.999727 

99.999727 
99.999727 


y,o(=) 


.000555497 

.000005638 

.000000546 

.000000391 

.000000010 

,000000012 

,000000012 


720 (m/s) 

250 

189.54 

198.90 

199.931 

200.009 

199.999576 

199.999612 

199.999612 


Round  E3 


2.5 

2.9987 

2.9999974 

2.9999997 

3.0000003 


y2QCm/s) 

250 

270.35 

298.31 

300.13 

300.02 

299.999242 

299.999324 

299.999324 


CjCl/m) 

.02 

-.0338 

-.0059 

.0082 

.00997 

.009998427 

.009998405 

.009998405 


c^Cs/m) 


.0 

.0096 
.0058 
.0122 
-.0024 
.000106786 
.000100384  j 
.000100384 


CP  ) 

n 

e(Pn)/e(Pn.i) 

std.  dev. 

205, 

.6972 

198.66 

. 00005 

.0050 

15,27 

.07688 

.0014 

3.33 

.21791 

.000645 

.20 

,05936 

.000157 

.0000013 

.00001 

.000000398 

,0000000023 

,00185 

.000000017 

.0000000023 

. 9999995 

.000000017 
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value  is  to  be  retained.  Since  it  is  the  user's  task  to  program  his 
particular  version  of  equation  set  (1)  oi  (2),  we  see  that  the  above 
feature  can  save  the  user  from  prograjming  many  versions  of  the  same 
equations,  the  versions  differing  only  in  the  nature  of  the  parair.eters 
involved.  If  the  version  programmed  contains  every  parameter  a 
reasonable  (or  only  slightly  unreasonable)  person  might  ever  want  to 
consider,  the  programmer  need  never  alter  his  program;  he  can  always 
suppress  unwanted  parameters  at  will. 

Of  course,  the  user  can  also  fix  any  paramic  at  a  nonzero  value. 
Consider,  for  example,  the  situation  where  some  of  the  input  paramic 
estimates  are  known  to  be  respectable,  ball-park  values,  while  the 
remaining  estimates  are  little  mere  than  wild  guesses.  There  is  no 
provision  in  FINLIE  for  weighting  the  paramic  estimates.  Thus  when 
the  data  are  especially  noisy,  FINLIE--in  Its  single-minded  effort  to 
decrease  e- -might  very  well  downgrade  an  excellent  estimate.  One 
way  to  avoid  (or  at  least  to  try  to  avoid)  this  difficulty  is  to  make 
two  computer  runs.  On  the  first  run,  all  highly  regarded  paramic 
estimates  are  held  fixed,  so  that  the  other  paramics  will  be  determined 
for  these  fixed  values.  The  fixed  and  determined  paramic  values  from 
this  first  run  then  ser^/e  as  the  estimates  for  a  second  run  in  which 
none  of  the  paramics  is  held  fixed. 

The  mechanics  of  informing  FINLIE  as  to  which,  if  any,  of  the 
paramics  is  to  be  held  constant  will  be  covered  later. 

In  Section  II,  we  discuss  in  more  detail  what  FINLIE  does  for 
the  user;  in  Section  III,  we  discuss  what  the  user  must  do  for  FINLIE. 


ir.  INSIDE  FINLIE:  WHAT  FINLIE  DOES  FOR  THE  USER 
We  rewrite  the  paramic  set  P  of  Eq.  (5)  in  the  form 

P  “  (Pjj  P2'  ’’’  (il) 

where  the  first  NR  x  N2  elements  of  P  are  the  initial  conditions  and 
the  remaining  N3  elements  ai'e  the  parameters. 

We  can  regard  P  as  a  point  in  an  N-dimensional  paramic  space  S. 
Then  e(P),  as  defined  by  Eq.  (7),  is  the  value  of  the  continuous 
scalar  point  function  e  at  point  P.  For  each  point  P  in  the  paramic 
space  S,  there  corresponds  a  single  value  eCP).  ^ FINLIE 's  task,  given 
a  starting  point  P^,  is  to  search  S  for  a  point  P  that  yields  a  minimum 

value  e(P).  (When  more  than  one  minimum  exists,  our  choice  of  starting 
point  Pp  usually  determines  whether  or  not  cCP)  is  the  desired 

absolute  minimum.) 
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The  fitting  process  carried  out  by  FINLIE  can  best  be  explained  in 
terms  of  a  single-round  situation.  Once  the  single-round  p/rocedure  has 
been  established,  it  will  then  be  relatively  easy  to  see  how  the 
process  can  be  extended  to  any  number  of  rounds. 

Hence  we  introduce  a  single-round  paramic  set  Q: 

Q  ^  ^2'  '*■ 

where 

N23  =  N2  +  N3,  (13) 

the  number  of  paramics  for  a  single  round.  Tlie  first  N2  elements  of 
Q  are  the  initial  conditions  and  the  remaining  N3  elements  are  the 
parameters.  For  our  sample  system  (3)  or  (4),  we  have  (for  any  one 
round) 


Q  ■  (ypo*  y2o'  “^i’  ‘^2^. 

Similarly,  we  introduce  a  single-round  version  of  e(P): 

N1 

VCO)  =  £  Z  "j„,  tyj„  -  yj(v«l’  f”’ 

m  j  =  l  ■>  J  J 

where  the  summation  on  m  is  over  the  measured  data  for  the  single  round, 
(For  Round  E2  of  Table  I,  for  example,  m  would  range  from  6  to  11,) 

Note  that  the  e  for  a  multi-round  situation^  Fq,  (7),  is  the  sum  of 
the  >'s  for  the  individual  rounds; 


e  = 


NR 


n=l 


(15) 


For  the  moment,  then--a  rather  long  moment,  lasting  until  Section 
TI  (G)--we  will  assume  that  FINLIE  is  handling  a  single-round 
situation:  only  one  set  of  initial  conditions  is  being  determined. 


A.  Condition  for  a  Minimum  y 


We  can  regard  Q  as  a  point  in  an  N73-dimensional  space  Sj.  A 

necessary  (though  insufficient)  condition  for  point  Q  to  yield  a 
minimum  value  of  y  is  that  the  gradient  of  y  at  that  point  be  the  zero 
vector: 


grad  y(Q) 


hM 


(16) 


Thus  FINLIE  must  seek  a  point  that  satisfies  all  N23  components  of  (16) 
simultaneously.  From  Eq.  (14),  we  see  that  at  any  point  Q 
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where 


=  -2  0^(Q)  [k-l,2,...N23] 


N1 


(17) 


m  jal 


and  where,  in  our  dimensional  notation, 

[«k>d  ■  (■'k'^’d 

[^jk^d  ~  ^d 

Thus  condition  (16)  can  be  vrritten  in  the  form 


0k  (Q)  =  0 


[k=l,2, . . .N23] 


The  N23  components  ^k  define  a  vector; 


t  (Q)  ■=  (0i(Q).  '  •  ®N23^Q^)si 


(19) 

(20) 
(21) 

(22) 

(23) 


which,  from  (16-17),  has  the  direction  of  the  negative  gradient  of  y 
at  point  Q;  that  is,  the  direction  in  which  the  rate  of  decrease  of  y 
is  greatest: 

t  =  -(1/2)  grad  y,  (24) 

^  is  a  vector  point  function  of  Q.  For  each  point  Q  in  the  paramic 
space  Sj^,  there  corresponds  a  unique  vector  Thus  FlNLIE's  search 
for  a  point  Q  that  yields  a  minimum  value  y(Q)  has  become  a  search  for 
a  point  Q  at  which  ^  is  zero. 

B ,  Influence  Coefficients 

The  partial  derivatives  D.,  in  (18)  are  sometimes  called  "influence" 

J 

(or  "sensitivity")  coefficients  because  they  reflect  the  influence  of 
the  paramics  on  the  solution  curves. 

To  satisfy  (22),  FINLIE  must  be  able  to  evaluate  the  influence 
coefficients  at  any  point  Q  for  each  independent  variable  value 
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The  manner  in 
equation  set. 


which  FINLIE  evaluates  D,.  (x  ,Q)  depends  on  which 

j  K  ni 

(1)  or  (2),  we  are  fitting  to  the  data. 


C.  Influence  Equations  for  System  (1) 


If  we  give  FINLIE  the  differential  equation  system  (1),  then  we 
must  also  give  FINLIE  a  system  of  differential  equations  for  the 
influence  coefficients.  Taking  the  partial  derivative  of  each  side  of 
(1)  with  respect  to  paramic  qj^,  we  have 


or,  assuming  that  the  order  of  differentiation  can  be  reversed. 


d  D., 

_ }JL 

dx 


n'=1.2,...N2  1 

[ k=l,2, .. .N23J 

The  system  (25)  is  subject  to  the  initial  conditions: 

D  (X  0)  =  i  ^ 

jk  lOotherwi! 


Lse 


(25) 


(26) 


(These  initial  conditions  merely  reflect  the  fact  that  the  influence 
coefficient  is,  by  our  definition,  8yj/3y^.p  and  hence  equals  one 


at  Xq.) 


The  paramics  affect  f,  (x,Y,C)  in  two  ways:  indirectly  through 

j 

their  effect  on  the  dependent  variable  vector  Y  and  directly  through 
the  parameter  vector  C.  Hence  (25)  can  be  rewritten  in  the  more 
cumbersome  but  (possibly)  more  revealing  form: 


[j=l,2,  ...  N2 
k=l,2,  ...  N23] 
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where 


subscript  C  indicates  that  x  and  vector  C  are  considered 
constant  in  taking  the  partial  derivatives  of  f^Cx.Y.C); 

subscript  Y  indicates  that  x  and  vector  Y  are  considered 
constant  in  taking  the  partial  derivatives  of  fjCx,Y,C). 

Thus,  by  "parsmic  differentiation”  we  obtain  an  auxiliary  system  of 
differential  equations  (27)  whose  solutions  are  the  influence  coeffi¬ 
cients  needed  to  fit  equation  set  (1) .  Note  from  (27)  that  those 
influence  equations  are  always  linear  in  the  influence  coefficients 
Djj^.  The  number  of  influence  equations  is 

NA  =  N2  X  N23  (28) 

The  user  must  include  his  version  of  system  (27)  in  the  FORTRAN 
subroutine  containing  his  version  of  system  (1). 

For  our  by-now-familiar  example,  system  (3),  we  have  NA  =  2  x  4=8. 
The  eight  influence  equations  for  system  (3)  are  shown  in  the  upper 
portion  of  Table  III,  where  (  )'  =  d(  )/dx. 

Recall  that  our  only  purpose  in  obtaining  the  influence  coefficients 
is  to  be  able  to  evaluate  6j^(Q),  Eq.  (18),  in  our  effort  to  satisfy 

condition  (22).  From  (18)  we  see  that  involves  only  for 

j=l  to  Nl;  that  is,  only  for  the  measured  variables.  Yet  Eqs.  (25-27) 

show  j  running  from  1  to  N2;  that  is,  over  all  the  dependent  variables. 

Do  we  have  more  influence  equations  here  than  we  need?  The  answer  is: 

no.  We  have  implicity  assumed  that  there  are  no  extraneous  dependent 

variables  in  system  (1) :  all  of  the  unmeasured  dependent  variables 

are  needed  to  solve  the  differential  equations  for  the  measured 

variables.  Hence  the  D.,  for  Nl  <  j  <  N2  are  also  needed. 

jk  ■> 

For  our  example,  system  (3)  with  Nl=l,  is  clearly  needed  to 
solve  the  differential  equation  for  y^.  Thus  each  is  also  needed, 
as  we  see  in  Table  III  (A).  (Cin  the  other  hand,  if  y^  had  been  the 
only  measured  variable  in  system  (3),  then  y^  would  be  an  extraneous 
variable  and  should  be  thrown  out.) 

The  mechanics  of  writing  and  submitting  the  influence  equations 
will  be  discussed  later.  FINLIE  will  automatically  assign  the  proper 
initial  conditions  (26)  and  integrate  the  influence  equations  simul¬ 
taneously  with  the  original  s>stem  (1)  to  obtain  yj (x^,Q)  and 

D., (x  ,Q)  at  each  x  . 
jk  m'^  m 
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Table  III.  Influence  Equations  for  System  (3)  and  for  System  (4) 
(A)  For  System  (3) : 

(DijV  E  3yj'  /ay^o  =  -^^2  ^*^21 

(0^2^'  -  /^y20  “  ^^22 

(Dij)'  =  3yj' /3Cj  =  -(y2  )D2j 

^^14^  ”  ^^1  ^^^2  ~  ~^^2  ^^24 

(D21)'  3y2'/3yio  =  -‘^i(1^2c2y2)D2j 

(D22)'  ^  ^^2' ^^^20  ""  ■‘^l ^^'*■^‘^2^2^ ^22 

(023)'  =  /:)c^  =  -0^(1+202/2^^25  ~ 

^^24^  “  ^^2  ^^^2  ~  “^1  ^^'*’^^2^2^ ^24  "  ^1^2 

where  (Djj)q  =  ^^22^0  ^  j^) ^  =  0  for  j  k 


Foi*  System  (4) ; 
“11  -  3/1/3/10  ' 
^12  3  3/1/3/20  ’ 

Di3  h  3yj/3Cj  = 
^14  ”  3y2/3C2 


=  -Cu-l)/(Cjy2Q  ) 

=  (b/c.  [1-u+c,  (x-xju] 


(U-1)/Cj  -  (x-Xq) 


(C)  Unneeded  Influence  Equations  for  System  C4) : 
“21  =  3/2/3/10  =  0 

^22  ”  3/2/3/20  “  ^/2'^/20^ 

2 

^23  “  "  ”^^2  ^ 

°24  "  3y2^^'^2  "  y2^ 


One  final  remark.  For  large  systems  with  many  paramics,  the 
exact  influence  equations  (27)  can  be  rather  cumbersome.  In  many 
cases,  certain  liberties  can  be  taken  with  the  influence  equations: 
expressions  can  be  approximated  by  simpler  ones,  the  effect  of  certain 
paramics  on  certain  terms  in  the  original  equations  can  be  ignored, 
etc.  If  done  with  care  and  judgment,  such  simpyf ications  will  have 
no  effect  on  the  final  answer:  the  same  p6int  Q  will  be  reached  with 
or  without  the  simplifications.  Note,  however,  that  discretion  is 
called  for.  If  the  user  has  any  doubts  as  to  the  merits  of  some 
modification  to  the  exact  influence  equations  (and  even  when  he  hasn't 
any  doubts),  his  safest  course  is  to  avoid  such  a  modification. 

D.  Influence  Equations  for  System  (2) 

If  we  give  FINLIE  the  solution  set  (2),  then  we  must  also  give 
FINLIE  the  influence  equations  obtained  by  differentiating  (2): 


3=1,2,  ...  N1 
k=k,2,  ...  N23 


(29) 


We  assume--as  with  system  (l)--that  there  are  no  extraneous  variables 
in  system  (2).  (For  (2),  this  means  that  the  initial  conditions  for 
all  of  the  unmeasured  dependent  variables  are  needed  to  evaluate  the 

expressions  for  the  measured  variables.)  However,  the  D.,  for 

J 

N1  <  j  <  N2  are  superfluous  and  should  be  ignored.  'IThus  the  number 
of  influence  equations  required  to  fit  system  (2)  is 

NB  =  N1  X  N23  (3C) 


To  fit  system  (4),  for  example,  (where  N1  =  1  and  N2  =  2),  the 

values  are  not  required  and  we  need  submit  only  four  influence  equations 

to  FINLIE.  These  equations  are  shown  in  Table  III  (B) .  FINLIE  will 

automatically  set  all  undefined  D, , 's  to  zero.  For  the  sake  of 

jk 

completeness,  expressions  for  the  unneeded  D,,  are  given  in  Table  III 

(C),  but  we  emphasize  that  these  latter  equations  should  not  be  given 
to  FINLIE.  Note  that  the  eight  expressions  for  in  Table  III  (B 

and  C)  do  indeed  satisfy  the  initial  conditions  indicated  in  part  A 
of  the  table. 


The  remarks  in  the  previous  section  on  the  possibility  of 
simplifying  the  influence  equations  apply  to  system  (29),  although 
here  the  urge  to  simplify  may  be  less  compelling. 
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E. 


An  Overview 


To  summarize  thus  far:  FINLIE  determines  the  values  of  y  ,  Q) 

1  ^ 

and  D.,  (x  ,  Q)  either 
jk  m  ^ 

(i)  by  numerically  integrating  a  system  of  N2  plus  NA  first- 
order  differential  equations  or 

(ii)  by  evaluating  a  system  of  N2  plus  NB  algebraic  or 
transcendental  expressions. 

Except  for  this  difference--but  what  a  difference  it  can  be  in  terms  of 
machine  time! --the  fitting  process  used  by  FINLIE  is  the  same  for  the 
two  equation  sets  (1)  and  (2). 

This  fitting  process  consists  of  adjusting  Q  until  it  satisfies 
condition  (22) .  Of  course,  it  would  be  pleasant  if  FINLIE  could  solve 
(22)  for  Q  in  some  direct,  one-step  fashion.  No  fooling  around  with 
Qq,  Qj^,  etc;  just  leap  in  and  solve  the  N23  equations  of  (22)  for  the 

N23  components  of  Q.  Unfortunately,  when  system  (22)  is  nonlinear  in 
one  or  more  of  the  paramics,  no  such  general  one-step  scheme  exists. 
Hence  FINLIE,  expecting  the  worst,  sets  out  to  solve  (22)  by  an 
iterative  process. 

Two  of  the  standard  iterative  techniques  are: 

(i)  differential  corrections  (alias  Taylor-series  lineariza¬ 
tion,  alias  Gauss  method,  alias  Gauss-Newton  method); 

(ii)  steepest  descent  (alias  gradient  search). 

FINLIE  uses  a  third  method,  due  to  Marquardt*,  which  is  a  blend  of  the 
first  two  methods,  retaining  the  best  features  of  each  and  avoiding 
their  disadvantages.  We  will  discuss  enough  of  the  differential 
corrections  and  steepest  descent  techniques  to  see  what  is  involved  in 
combining  the  two. 


F.  Differential  Corrections  in  Space 


For  each  point  Q  in  there  corresponds  a  position  vector 
Let  be  the  vector  from  point  Q  to  point  Q: 


I 


=  5  .  (Iq,.  •  •  ■  ''■lN23ls,  ^ 


(-M) 


*  r<ec  the  Bibli  ogvcxphij ,  Part  A. 


. . .  mnijrswfftwwwBW 


In  the  differential  corrections  technique,  we  approximate  the  basic 
condition  (22)  by  a  system  of  equations  (to  be  derived  in  the  next 
paragraph)  that  is  linear  in  the  incrranents  Aq^^.  We  can't  solve  (22) 

for  Q,  but  given  a  point,  say  Qq,  we  can  solve  the  approximate  condi¬ 
tions  for  an  approximate  increment  vector  This  increment  is  then 

added  to  to  reach  the  next  way- station  on  our  trek  to  Q: 

$j  •  5o  *  iSo 


Point  Qj  is  ^.n  improvement  over  point  Qq  if  yCQj)  is  less  than  yCQq)- 

But  improvement  or  not,  the  differential  corrections  method  plows 
ahead,  u"ing  Qj^  to  re- solve  the  approximate  equations  for  a  new  incre¬ 
ment  i’h®  process  continues  in  this  manner  through  a  series  of 

points  until  a  specified  convergence  criterion  has  been  met  or  a 
specified  number  of  iterations  have  been  performed  or  some  numerical 
catastrophe  occurs. 

The  desired  approximation  to  condition  (22)  can  be  obtained  by 
expanding  y.  and  D.,  in  Taylor  series  about  point  Q.  We  have 

1  J  K 


y-  (x  ,6)  =  y.  (x  ,Q)  + 


N23 

D.  (X , 
^  jn  m' 


Q)  • 


+  (higher-order  terms) 


D.,  (x  ,Q)  =  D.,  (x  ,Q)  +  (higher-order  terms) 
jk  m’^  j  k  m'^ 


(33) 


(34) 


We  assume--an  assumption  that  is  not  always  valid--that  Q  is  close 
enough  to  Q  to  permit  us  to  ignore  the  higher-order  terms  in  (33)  and 
(34).  Then  from  definition  (18),  we  have 


N1 


r 


m  j  =  i 


w. 

jm 


N23  "I 

-  y-  (x  ,Q)  -  y  .  D.  (x  ,Q)  •  Aq  •  D.,  (x  , 
jm  ■'j  jn  m’^  ^n  j  jk  m’ 


Q) 


6^(Q)  - 


Tv 


23 


By  rearranging  the  sums,  we  obtain 


e^CQ)  ^  e^(Q)  - 


(35) 


where 


“kn")’ 


V  ^  j"*  J 


,  (x  ,Q)  •  D.  (x  ,Q) 
k  m'^  jn  m'^ 


(36) 


Thus  the  conditions  B]^(Q)  ~  0»  which  hold  at  a  point  Q  where  y  is  at  a 
minimum,  are  replaced  by  the  conditions 

I 

(37) 

_J 

[k  =  1,  2,  ...  N23] 

which  are  applicable  to  points  in  the  vicinity  of  Q. 

The  quantities  Uj^j^(Q)  have  at  least  four  interesting  properties: 

K„ld  ” 


a  ,  =  a, 
nk  kn 


“kk  ^  ® 


(38) 


a  n, ,  >  a  , 
nn  kk  nk 


The  first  three  properties  follow  at  once  from  definition  (36);  the 
fourth  is  a  consequence  of  Holder's  Inequality  (alias  Cauchy's,  alias 
Schwarz's,  alias  Buniakovski ' s  Inequality).  In  general,  we  regard 

as  the  (k,n)-th  element  of  an  N23  by  N23  s>Tnmetric  matrix  a. 

In  matrix  form,  (37)  becomes 


kn 


[a(Qj  •  =  e'^(Q)] 


1 


(39) 
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where  the  superscript  T  (for  Transpose)  denotes  a  column  vector  and  the 
subscript  Sj  indicates  that  all  components  are  in  the  N23-dimensional 

space  Sj.  For  either  of  our  examples,  system  (5)  or  (4),  (39)  becomes 


/“ll 

“l2 

“l3 

“21 

“22 

“23 

“31 

“32 

“33 

\“41 

“42 

“43 

‘24 


‘34 


“44/ 


) 

r\ 

^^20 

A"l 

Vc2  / 

V./ 

(40) 


System  (39)  is  linear  in  the  increments  Aqj^;  hence  the  process  of 

solving  for  these  increments  is  routine  work  for  the  computer.  (We 
assume  that  a  solution  does  exist;  this  amounts  to  assuming  that  the 
determinant  of  matrix  a  is  not  zero.) 


■in«  corrections  process,  then,  consists  of  substitut¬ 
ing  Qq  in  (39),  solving  for  substituting  in  (39)  the  point  Q 

obtained  by  the  vector  addition  solving  for  AQj,  elc. 

Unfortunately,  even  when  this  process  converges  to  some  point 

absolute  minimum 

Y.  Condition  (2^)--which  is  approximated  by  the  matrix  equation  (39)-- 
guarantees  only  that  its  solution  point  Q  will  yield  some  Relative 
extremum  value  of  y.  Space  Sj  could  be  teeming  with  points  of  local 

iTrT^rt  J^ese  extremum  points,  including  the  one  we  seek, 

a  sort  of  black  hole  in  space  S^,  capable  of  drawing  a  nearby  search 

party  into  its  core.  The  particular  black  hole  into  which  we  are  drawn 
depends  mainly  on  where  we  start  in  space  . 

Differential  Corrections  in  Space  S 


So  far  in  Section 
For  this  situation,  the 
matrix  equation  (39). 


II,  we  have  assumed  single-round  data,  NR  = 
differential  corrections  technique  led  to 


1  . 


Ei  T  situation  of  Table  1.  For  each  round 

ti  11  -  1.  2,  3),  FINLIE  forms  a  vector  and  a  matrix  by  Eqs.  (18) 

and  (36)  respectively,  using  the  Q  and  m  indicated  below: 


Range  of  Subscript  m 

Round  _ Point  Q _ in  Eqs.  (x8)  and  C36j 


El 

^^lO^El* 

^y20^El’ 

*^2 

1 

to 

5 

E2 

^^10^52’ 

o 

m 

ISJ 

^1* 

6 

to 

11 

E3 

^^lO^ES’ 

'•y20^E3* 

"r 

"'2 

12 

to 

16 

In  the  8-dimensional  space  S  associated  with  the  eight  paramics  p^  of 
Eq.  (10),  the  vectors  take  the  form; 


t 

t 


El 

E2 

E3 


[(BiIep  (62)ei.  0,  0,  0,  0,  (B^)^^, 

[0,  0,  (ei)g2»  f^2^E2'  ^^3^E2' 

[0,  0,  0,  0,  CB2)e3’  ^S^E3* 


^‘^4^E2^S 


‘■^4^E3^S 


(41) 


Similarly,  in  space  S  the  matrix  a  for  round  El  expands  to: 

A^iPeI  *-°‘12^E1  0  0  0  0  (otj3)pj  f“i4)r.;i\ 


El 


(ci2i)pi  (a22^EI 


0 

0 

0 

0 


0  0  0  0  (a,^),,. 


0 

0 

0 

0 


0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 


*23' 

0 

0 

0 

0 


(a  _ , )  „  ^  (a_  *) , 


0000 


*■“24^  El 


0 

0 

0 

0 


^°‘34^ni 


31' El  "3 2' El 

V°‘4pEl  ^'^42^E1  0  0  0  0  (o<43)j;j  f^44JE 


(42) 


El/ 


with  similar  expressions  for  Up,,  and 


Since  the  multi -round  e  to  be  minimized  is  the  sum  of  the  single- 
round  >'s,  EINLIE  obtains  the  multi-round  version  of  E.q.  (40)  by  slimming 
--in  space  S--the  three  single-round  .  vectors  of  Eq.  (41); 


(43) 


t  =  [tgj  +  '?g2  ^e3^S 


and  the  three  single-round  matrices: 

A  =  [ogj  +  ag2  + 

The  desired  multi-round  matrix  equation  is  then 

[A(P)  • 

A  detailed  form  of  this  equation  for  our  three-round  sample  system  is 
given  in  Table  IV;  the  generalization  to  any  number  of  rounds  can  be 
easily  visualized. 

The  N  by  N  symmetric  matrix  A  will  always  contain 

NR  X  (NR-1)  X  N2  X  N2 

zeroes  distributed  among  the  off-diagonal  elements  of  all  but  the  last 
N3  rows  and  columns.  Let  be  the  (k,n)-th  element  of  matrix  A.  As 

in  Eqs.  (38),  we  have 

^^kn^d  ~  f^Pk^n^  ^d 

^nk  ~  ^kn  1 

(46) 

\k '  “ 

a  a, ,  >  a  ,  ^ 
nn  kk  nk 

Similarly,  if  bj^  denotes  the  k-th  component  of  vector  then 

[b,]d  -  [P,-‘ld  07) 

We  have  taken  some  pains  to  distinguish  between  the  multi-round 
paramics  p^^  and  the  single-round  paramics  which  for  our  three-round 

sample  systems  take  the  form 

P  =  0*^20^ El'  ^^10^20^E2'  *^^10^20^E3’  ‘^l’  ^2^ 

Q  =  ^yio’yzO”  ‘^l’  ^2^ 
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able  IV.  Matrix  Equation  for  System  (3)  or  (4),  Given  Data  from  Three  Rounds 


The  chief  reason  for  taking  these  pains  is  that  the  FINLIE  user  must 
himself  make  this  distinction  in  a  multi-round  situation.  The  FINLIE 
input  arguments  (to  be  discussed  later)  are  defined  in  terms  of  the  N 
paramics  pj,,  but  the  influence  equations  submitted  to  FINLIE  must  alway 

be  written  in  tenus  of  the  N23  (=  N2  +  N3)  paramics  qj^.  The  values  of 

the  initial  conditions  may  change  with  the  round,  but  the  influence 
equations  themselves,  like  the  original  equations  (1)  or  (2)  on  which 
they  are  based,  remain  the  same.  Thus,  regardless  of  the  number  of 
rounds,  there  will  always  be  N23  influence  coefficients,  defined  in 
terms  of  the  N23  paramics  qj^,  and  there  will  always  be  NA  (=  N2  x  N23) 

or  NB  (=  N1  X  N23)  influence  equations  (depending  on  whether  the  user 
is  working  with  system  (1)  or  system  (2)), 

H.  Differential  Corrections  in  Space  S 

If  the  paramics  pj^  are  not  all  of  the  same  dimension,  our  paramic 

space  S  is  a  hodgepodge;  a  salmagundi,  a  gallimaufry,  an  olla-podrida 
of  units.  Certain  computational  advantages  can  be  obtained  by  working 
in  a  space  S  in  which  the  paramics--and  hence  the  components  of  grad 
e--are  nondimensional .  (The  advantages  of  S  are  especially  compelling 
in  the  steepest  descent  technique,  some  of  whose  properties  are  not 
scale- invariant . ) 

To  achieve  the  desired  paramic  transformation  from  S  to  S,  we  note 
from  Eqs.  (46)  that 


or 


Pv]  =  1 


(48) 


That  is,  the  bracketed  quantity  in  (48)  is  nondimensional.  Thus  the 
paramic  transformation 


Pk  =  (a^ 


kk-'  Pk 


(49) 


creates  the  desired*  paramic  space  S.  The  elements  of  A  and  B  in  S  are 

*Fr>om  Eq.  (47),  we  see  that  the  product  bj^pj^  is  also  nondimensional. 
Thus,  the  transformation  ^ 

Pk  "  ^kPk 

seems  appealiruj;  it  would  lead  to  a  space  in  which  all  components  of  ? 
are  unity,  'flic  .qrpeal,  howei’er,  is  illusory,  Tt  would  not  be  vcrnj 
wise  to  use  as  scale  factors  the  very  quantities  h.  that  we  arc  tvyi.jij 


to  drive  to  zero, 
never  zero  (sec  Ea .  (46)) 


The  scale  factors 


on  the  other  hand,  are 
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(SO) 


“jk  =  “kk>''‘»jk 


These  space  S  components  have  the  following  admirable  features; 


(i) 

(ii) 


p  ,  a.,  and  b,  are  nondimensional ; 

K  J  K  K 

the  diagonal  elements  of  matrix  A  are  unity; 


\k  =  1 


(52) 


(iii)  the  off-diagonal  elements  of  A  satisfy  the  inequality 

(53) 


-1  <  a.,  <  1 

Finally,  the  form  of  matrix  equation  (45)  is  unchanged: 


[  A(P)  •  =  t  ^{P)  ] 


(54) 


the  subscript  S  serving  to  remind  us  that  all  components  are  now  in 
the  scaled  naramic  space.  Ninety-nine  percent  of  the  labor  in  solving 
(54)  for  AP  is  usually  expended  in  inverting  matrix  A.  Use  of  the 

scaled  components  a^j^  tends  to  increase  the  accuracy  of  the  matrix 

inversion  process. 


Note  that  each  scale  factor  (Sj^j^)  ^  in  (49)  is  a  function  of  point 

P,  the  curi'ent  set  of  paramic  values.  Hence  each  time  the  paramics 
are  up-dated,  a  new  transformation  must  be  made:  a  new  S  space 
created.  "his  i  ^no  big  problem  for  a  computer.  FINLIE  nandles  the 
scaling  to  space  S  and  back  again  to  the  user's  space  S;  the  process 
is  automatic  and  invisible  (in  computer  jargon,  "transparent")  to  the 
user . 


I .  Steepest  Descent 

Consider  a  given  point  P^  and  the  corresponding  vector  B(Pp) 

proceeding  from  that  point.  Recall  that  5  at  any  point  is  a  vector 
in  the  direction  of  the  negative  gradient  of  c  at  that  point.  Hence, 

provided  that  the  magnitude  of  t(P„)  is  not  zero  (if  it  were,  P  would 

,  o  o 

be  the  desired  solution  P),  ^(Pq)  is  the  steepest  descent  vector  for 
point  P^;  a  vector  in  whose  direction  c(P)  will  decrease  most  rapidly 
(at  least  at  first)  as  we  move  away  from  P^.  Let  Pj  be  any  other  point 
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in  this  steepest  descent  direction: 

^  ^  '  t(Fo)  (55) 

where  h  is  a  nondimensional  positive  constant. 

There  always  exists  a  range  of  h  values,  0  <  h  <  which 

the  point  obtained  by  (55)  is  an  improvement:  ^(Pj)  e(Po)-  The 

steepest  descent  method  determines  the  optimum  h  in  this  range;  the 
value  of  h  for  which  e  is  a  local  minimum  along  the  vector  ^(Pq) .  This 
can  be  done  by  evaluating  Pj  and  e(Pj)  for  a  series  of  h  values: 

ho  <  hj  <  h,  <  .  .  .  Presumably,  for  a  while  e  will  decrease  with 

increasing  h.  As  soon  as  an  h  is  found  for  which  the  e  has  increased, 
the  (approximately)  optimum  h  for  point  Pp  can  be  determined  by 
interpolation. 

Given  the  new  point  Pj  based  on  this  optimum  h,  the  next  point 
P2  will  lie  in  the  direction  of  steepest  descent  from  P^;  that  is, 
along  the  new  vector  ^(P^) .  Another  optimum  h  must  be  determined  to 
obtain  P2-  And  so  on  to  P. 

The  diff iculty^with  this  approach  is  that  in  the  neighborhood  of 
the  solution  poipt  P,  where  |?j  is  nearly  zero  and  yet  we  are  not  quite 
close  enough  to  P  to  be  able  to  quit  with  honor,  further  progress 
is  painfully  slow.  Often  the  sampling  size  on  h,  the  Ah  intervals, 
must  be  shortened  beyond  all  endurance  in  an  effort  to  find  a  Pj 
for  which  e(Pp  <  e(PQ).  Ingenious  variations  on  the  basic  steepest 

descent  theme  have  lessened  but  not  removed  this  difficulty. 


J.  Marquardt  Interpolation  in  Space  S 

The  two  fitting  techniques  we  have  discussed  so  far  are: 

(i)  differential  corrections,  which  in  space  S  is  based  on  matrix 
equation  (54) ;  this  equation  has  the  component  form 


N 

n=l 

[k=l,2,. .N1 


(ii)  steepest  descent,  based  on  the  vector  equation  (55),  which 
in  space  S  nas  the  component  form 
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hb^CP) 


(57) 


^Pk 


[k=1.2....N] 


Comparing  these  two  methods,  we  note  that; 


(a)  Far  from  the  solution  point,  the  steepest  descent 
technique  is  superior.  It  must  proceed  so  as  to  decrease 
e,  whereas  the  differential  corrections  method  is  under 
no  such  compulsion  and  is  likely  to  lead  us  into  strange 
pastures . 

(b)  Close  to  the  solution  point,  the  differential  corrections 
method  is  superior.  It  converges  rapidly  in  the  very 
region  where  the  steepest  descent  technique  languishes. 

Marquardt*  has  proposed  an  interpolation  between  the  two  methods; 
a  technique  that  behaves  like  the  steepest  descent  when  we  are  far  from 
the  solution  and  like  the  differential  corrections  method  when  we 
enter  a  neighborhood  in  which  the  higher-order  terms  in  Eqs.  (33)  and 
(34)  are  negligible. 


To  achieve  this  interpolation,  a  positive  nondimens ional  constant 
A  is  added  to  each  diagonal  element  of  the  scaled  matrix  A.  That  is, 
the  system  (56)  is  replaced  by 


(58) 


where  ~ 

a,  =  1  +  A  when  k=n 

kn 

=  a,  when  k/n 
kn 


(59) 


System  (58)  is  the  bedrock  upon  which  the  FINLIE  fitting  process  rests. 
Note  the  behavior  of  this  system  as  a  function  of  A; 


(a)  As  A-K),  system  (58)  clearly  reverts  to  the  differential 
corr^tions  system  (56). 

(b)  As  A->  the  diagonal  terms  of  system  (58)  dominate  and 

the  system  degenerates  into  N  uncoupled  equations  of  the  form 

(1  +  A)  \ 


*Sec  the  Bibliography,  part  A, 
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or,  since  by  assumption  X>>1, 

Apj^  =  \ 

Comparing  (60)  with  (57),  we  see  that  for  large  A^values,  system  (58) 
simulates  the  steepest  descent  approach  with  h=X”  .  That  is,  for 
X>>1,  (58)  will  take  us  to  a  new  point  a  rather  short  distance  from  the 
current  point  P  in  the  direction  of  the  negative  gradient. 

Marquardt  has  suggested  an  algorithm  for  determining  a  suitable 
value  of  X  for  each  iteration;  that  is,  for  each  step  to  Pj,  Pj  to 

?2t  etc.,  on  the  path  to  the  desired  solution  point  P.  This  algorithm 

(with  a  few  very  minor  "refinements")  has  been  incorporated  into 
FINLIE.  The  basic  scheme  is  as  follows. 

For  the  first  iteration,  P.  to  P  ,  FINLIE  assigns  a  tentative  value 
to  X:  u  i 


(starting  X)p^  =  Xj^^  =  0.001  (61) 

Let  Pj^^  denote  a  candidate  for  point  Pj,  obtained  by  solving  (58)  with 
P.Pq  and 

’  ^0  *  ‘■“('’o-hA' 

The  basic  test  that  any  point  P  should  pass  is  that  it  be  an  improve¬ 
ment  over  the  current  point: 


e(P)  <  e(PQ)  (63) 

If  Pjy^  satisfies  test  (63),  then  FINLIE  returns  that  point  to  the  user 
as  the  updated  point  P,  and  is  ready  to  start  the  next  iteration,  P,  to 

P  *  i 

If  P-^  fails  test  (63),  then  FINLIE  must  take  a  smaller  step  in  a 

more  propitious  direction.  This  can  be  accomplished  by  increasing  X. 
That  is,  FINLIE  re-solves  system  (58)  with  P=P„  as  before,  hut  with 
X  increased  to,  say,  — — — - 


X 


IB 


10  X 


lA 


(64) 


(Note  that  in  re-solving  the  system  (58),  the  elements  (kji^n)  and 

do  not  have  to  be  re-evaluated.  They  depend  only  on  the  current 

point  and  thus  are  evaluated  only  once  each  iteration.)  The  new 
increment  vector  for  X^^g  yields  the  new  candidate  point: 
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(65) 


If  Pjg  satisfies  test  (63),  then  FINLIE  returns  this  point  to  the 

user;  if  fails  test  (63),  then  FINLIE  increases  X  again  by  a  factor 

of  ten,  and  so  on.  Sooner  or  later,  an  acceptable  candidate  will  be 
found : 


H-  Ip  (Pq.io^^ia'J  (66) 

where  n  is  zero  or  a  positive  integer. 

The  cost-conscious  reader  may  ask:  if  P^^  fails  test  (63),  why  not 

skip  over  a  possibly  long  line  of  rejected  candidates  by  increasing  A 
by  some  factor  much  larger  than  ten?  This  should  get  us  to  an  accept¬ 
able  candidate  point  at  once  or  at  least  in  fewer  trials.  True,  but 
the  general  principle  is  this:  the  larger  the  A,  the  smaller  the 
progress  we  are  making.  Hence  we  don't  want  FINLIE  to  use  a  A  "very 
much"  larger  than  needed  to  satisfy  test  (63).  It  is  not  worth  the 
effort  to  find  the  optimum  A  for  each  iteration,  but  by  testing  after 
each  ten-fold  increase  in  A,  FINLIE  will  not  grossly  exceed  that  op¬ 
timum.  (Indeed,  a  case  could  be  made  out  for  merely  doubling  A  each 
time  an  increase  is  required.) 

The  only  way  in  which  the  second  and  subsequent  iterations  differ 
from  the  first  is  in  the  formula  FINLIE  uses  for  determining  the  start¬ 
ing  A  value  for  the  iteration: 

(starting  A)p  to  P  “  0.1  x  (final  A  value  used  to  (67) 
n-1  n  produce  point  P^  ^  in 

the  previous  iteration) 

That  is,  FINLIE  always  decreases  the  current  value  of  A  by  a  factor  of 
Len  at  the  start  of  each  new  iteration.  This  decrease  is  an  essential 
part  of  the  A  manipulation.  When  all  is  going  well,  FINLIE  will  have 
no  need  to  increase  A;  thus  rule  (67)  will  insure  that  A  goes  to  zero  - 
and  hence  that  the  process  approaches  the  differential  corrections 
technique  -  as  FINLIE  approaches  the  solution  point  P. 

A  typical  set  of  A  values  encountered  in  the  course  of  some 
hypothetical  fit  (not  our  familiar  examples,  (3)  and  (4))  is  shown  in 
Table  V.  The  reader  can  infer  from  these  A  values  the  fleeting 
existence  of  rejected  candidate  points.  Thus,  to  get  from  P^  to  P^, 

FINLIE  clearly  had  to  solve  system  (58)  six  times:  for  A  =  10“'’  (that  is, 

one-tenth  the  previous  A),  10‘^,  lO’^,  10'\  10°  and  10^  (the  A  value 
that  produced  a  successful  candidate).  Similarly,  to  get  from  P^  to  P^, 
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TABLE  V.  Typical  X  Values  During  a  Fit 

Iteration 

X  value  returned 
by  F INLIE  at  the 
end  of  the  iteration 

No.  of  times  system  (58) 
must  have  been  solved  by 
FINLIE 

Pq  to  Pj 

10“^ 

Pi  P2 

10“^ 

1 

^2  P3 

10^ 

6  (for  X-IO'"^,  10’^....  10^] 

P3  to  P^ 

10° 

1 

P4  to  P3 

10’^ 

1 

P5  ^6 

10'^ 

2  (for  X=10'^,  IQ'^) 

P6  P7 

0 

1 

K> 

1 

P?  to  Pg 

10"^ 

1 

Pg  to  Pg 

lO'"* 

1 

‘’g  “  '’10 

0 

1 

</l 

1 

PlO  to  Pjj 

10'^ 

1 
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FINLIE  must  have  solved  (S8)  twice:  for  X»=10”^  and  10  .  Thereafter, 

the  fitting  process  seemed  to  get  back  on  the  track  and  X  decreased 
steadily.  Without  Marquardt’s  X  in  the  system,  it  is  likely  that  the 
search  represented  by  Table  V  would  have  gone  astray  after  point  P2  and 
come  to  some  abrupt  and  ignoble  conclusion. 


K.  Convergence  Cr.iterion 

The  question  arises:  wnen  can  the  user  accept  a  point  returned  by 
FINLIE  as  being  "close  enough"  to  the  desired  solution?  One  possible 
answer  is:  when  FINLIE  tells  him  he  can.  At  the  end  of  each  iteration, 
FINLIE  returns  to  the  user  a  flag  whose  value  indicates  whether  or  not 
the  returned  point  has  satisfied  a  built-in  convergence  criterion. 

(This  flag  will  be  discussed  in  section  111(C).) 

The  convergence  criterion  installed  in  FINLIE  is  as  follows.  Let 
P  j  and  P  be  any  two  consecutive  points  returned  by  FINLIE:  the 
efia^  ■points^of  two  consecutive  iterations.  Then  FINLIE  will  signal 
convergence  at  point  P^  if  and  only  if 

0.99999  <  ^ 

The  right-hand  portion  of  this  double  inequality  is  essentially 
inequality  (63)  and  hence  is  always  satisfied,  thanks  to  the  Marquardt 
X  feature.  The  left-hand  inequality  in  (68),  however,  constitutes 
an  arbitrary  definition  of  convergence:  namely,  that  the  percent 
change  in  e  has  dropped  below  0.001. 

As  an  example  of  criterion  (68)  in  action,  consider  the  search 

summarized  by  Table  II.  The  values  of  CRHe(P  )/e(P  ,)  listed  in  the 

'  n  n-1 

next  to  last  column  of  that  table  jump  about  erratically  (always 

between  0  and  1,  of  course)  before  the  criterion  is  satisfied  at  point 

P^.  The  sudden  transition  from  the  value  of  CR  at  P,  to  its  value  at 

Py  is  not  typical.  In  searches  based  on  more  realistically  inaccurate 

measured  data,  CR  will  often  be  close  to  -  and  monotonically  approach  - 
the  value  1  over  the  final  few  iterations. 

Note  that  (68)  is  only  a  measure  of  convergence  to  a  local 
minimum.  We  have  said  it  before,  but  it  bears  repeating:  there  is  no 
guarantee  that  the  point  P  satisfying  (68)  will  yield  the  desired 
absolute  minimum  e. 

Of  course,  the  user  need  not  accept  definition  (68);  he  can  ignore 
the  FINLIE  convergence  flag  and  impose  his  own  convergence  test  on  the 
data  returned  by  FINLIE  after  each  iteration. 
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L.  Estimated  Errors 


In  addition  to  computing  the  estimated  standard  deviation  of  the  fit 

1  h 


s  = 


(69) 


FINLIE  computes  s^,  the  estimated  standard  deviation  of  paramlc  pj^, 
k=l, 2, . . .N. 


For  linear  least- squares,  the  conventional  formula  is 

"Ic  ■ 

where 


(70) 


A, ,  =  the  k-th  diagonal  element  of  the  inverse 
^  of  the  unsealed  matrix  A. 

Note  that  while  s  is  nondimens ional,  Sj^  has  the  same  dimensions  as 

Pk* 

For  nonlinear  least-square  fits,  Eq.  (70)  should  be  viewed  with  a 
healthy  suspicion.  Indeed,  CelmiipS  (Ref.  33  in  the  Bibliography)  points 
out  that  even  in  the  linear  case,  the  equation  should  be  applied  only 
in  "very  limited  special  cases."  Unfortunately,  the  alternative  formula 
that  he  develops  for  Sj^  is  a  rather  complicated  one  involving  second- 

order  derivative  terms  -  terras  that  so  far  we  have  managed  to  avoid. 

The  inclusion  of  these  terms  would  mean  more  work  not  only  for  FINLIE  - 
which  would  be  acceptable  -  but  for  the  user,  who  would  have  to  derive 
and  program  some  possibly  horrendous  expressions.  The  labor  here 
seems  out  of  proportion  to  its  reward,  since  the  "crude"  error 
estimates  provided  by  (70)  are  usually  not  all  that  crude  when  the 
search  has  converged  to  the  proper  point.  Hence  FINLIE  returns  these 
estimates  to  the  user  and  the  user  is  expected  to  provide  his  own  grain 
of  salt. 

(Note  that  (70)  uses  only  the  diagonal  elements  of  the  inverse 
matrix.  In  some  situations,  all  of  the  elements  of  A"  are  useful 

-12 

for  error  analysis.  In  these  special  situations,  A  s  can  be  regarded 
as  the  variance- covariance  matrix.  However,  for  nonlinear  least 
squares,  we  are  pushing  our  luck  in  making  use  of  the  diagonal  elements; 
to  try  to  assign  any  significance  to  the  off-diagonal  elements  would 
really  be  folly.) 

Recall  that  FINLIE  transforms  the  elements  of  matrix  A  to  the 
scaled  space  S,  Eq.  (50),  and  then  replaces  the  diagonal  elements  by 
1+X.  Hence  FINLIE  actually  obtains  the  paramic  error  estimates  by  the 
relation 
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where 


s 


k 


®kk 


is 


s 


(71) 


the  k-th  diagono^l  element  of 

the  inverse  of  matrix  Eq.  (58) 


(I  felt  there  should  be  some  compensation  in  the  error  estimate 
formula  for  the  presence  of  Marquardt's  X  in  the  equations.  By  a  chain 
of  nonrigorous  reasoning,  I  was  thus  led  to  insert  the  (1+X)  factor  in 
(71).  Since  X«1  for  a  good  fit,  (1+X)  seems  relatively  harmless 
sitting  there.) 


M.  The  Composition  of  FINLIE 

So  far,  the  word  FINLIE  has  denoted  an  apparently  monolithic 
program.  Actually,  for  reasons  that  seemed  persuasive  at  the  time, 
FINLIE  was  written  as  an  assemblage  of  six  linked  FORTRAN  subroutines: 

DUBLIN,  LONDON,  PARIS,  BONN,  MATINV,  MERSO 

only  one  of  which  -  DUBLIN  -  is  called  by  the  user.  "FINLIE",  then, 
is  merely  a  convenient  name  for  an  ensemble  of  six  subroutines. 

[FINLIE  is  also  the  name  of  a  permanent  file  (in  Update  format) 
stored  on  the  front  end  of  BRL's  Control  Data  Corporation  computer 
system.  (At  BRL,  this  system  consists  of  two  linked  mainframes:  the 
CYDER  170/Model  173  and  the  CYBER  70/Model  76.)  File  FINLIE  contains 
five  of  the  six  subroutines:  all  but  MATINV,  which  is  already  available 
from  a  systan  library.] 

The  relationship  between 

(i)  the  user's  program  that  calls  FINLIE, 

(ii)  FINLIE 

and  (iii)  the  user's  subroutine  defining  his  equations, 

and  the  inter-relatiorship  of  the  six  subroutines  that  constitute 
FINLIE  are  all  indicated  in  Figure  1.  A  vertical  bar  between  two 
subroutines  in  the  figure  indicates  that  the  upper  subroutine  calls 
the  lower  one. 
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All  six  subroutines  of  FINLIE  are  listed  in  the  Appendix.  Only 
four  of  the  six  -  the  four  "cities"  -  were  written  by  the  author;  the 
other  two  Cnamely,  MATINV  and  MERSO)  are  general-purpose  subroutines 
to  be  discussed  shortly.*  With  some  minor  exceptions  in  subroutine 
MERSO  (these  will  be  spelled  out),  the  FORTRAN  used  in  FINLIE  is  a 
"more  or  less  standard"  version  of  FORTRAN  IV  (alias  FORTRAN  4,  alias 
FORTRAN  66). 

Converting  FINLIE  to  a  later  model  FORTRAN  -  say,  FORTRAN  77  - 
should  be  relatively  uneventful.  One  possible  difficulty  is  as  follows. 
FINLIE  was  written  for  a  compiler  that  automatically  retains  the  values 
of  entities  defined  within  a  subroutine  but  not  linked  to  the  calling 
program.  For  such  a  compiler,  subsequent  calls  to  the  subroutine  will 
find  the  previous  values  waiting.  However,  in  FORTRAN  77  the  SAVE 
statement  is  available  for  specifying  what  if  anything  is  to  be  retained; 
hence  some  FORTRAN  77  compilers  may  not  automatically  retain  local 
values.  In  that  case,  it  may  be  necessary  to  SAVE  the  arrays  ALPHA  and 
BETA  in  subroutine  DUBLIN. 

DUBLIN  is  the  interface  between  the  user  and  FINLIE.  The  user  must 
write  the  FORTRAN  program  that  calls  subroutine  DUBLIN  with  the 
required  input  data.  Hereafter,  we  will  refer  to  (and  think  of)  the 
user's  calling  program  as  a  main  program,  although  it  could  itself  be  a 
subprogram.  Each  time  that  DUBLIN  is  called  by  this  main  program, 

DUBLIN  activates  the  other  subroutines  of  FINLIE,  causing  one  iteration 
of  the  search  procedure  to  be  carried  out.  That  is,  if  the  main  program 
submits  point  to  DUBLIN,  DUBLIN  will  return  to  the  main  program 

the  next  point  P^^.  Information  and  advice  on  writing  the  main  program 

and  in  particular  on  calling  DUBLIN  will  be  given  in  overwhelming  detail 
in  Part  III. 

Subroutines  LONDON,  PARIS  and  BONN  are  buried  within  FINLIE,  so  that 
their  individual  purposes  should  be  of  little  significance  to  the  user. 
However,  the  following  features  of  PARIS  can  be  noted  from  Fig.  1.  If 
the  user  is  fitting  a  set  of  differential  equations,  PARIS  calls  a 
numerical  iiitegration  subroutine  MERSO  (of  which  more  will  be  said 
shortly)  and  MERSO  in  turn  calls  the  subroutine  -  written  by  the  user 
and  arbitrarily  labelled  ROME  in  the  figure  -  that  defines  the 
differential  equations  to  be  fitted.  On  the  other  hand,  if  the  user 
is  fitting  algebraic  or  transcendental  equations,  PARIS  calls  the  user's 
equation-defining  subroutine  ROMA  directly.  Both  ROME  and  ROMA  must 
get  additional  information  from  PARIS  through  the  labelled  COMMOfT 

reviewer  of  thin  piz/'t’r  ijuestioned  the  imf^lioation  that  DUBLIN, 

LONDON,  PARTS,  atiJ  BONN  ai’e  the  only  "oit  iar'  in  the  oextup  let  of 
auhroutines.  He  went  so  far  as  to  consult  an  atlas  to  see  if  there  is  a 
tomx,  a  village,  a  hamlet  or  a  crossroads  by  th'  n^ane  of  MATINV  or 
MFRSO  somewhere  i?i  the  Appiia^ently  there  isn't. 
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block  NAPLES.  ROME  and  ROMA  may  require  additional  information  from  the 
user's  main  prograjn,  this  can  be  passed  through  the  labelled  COMMON 
block  CAIRO.  Abundant  details  on  writing  ROME  and  ROMA  and  on  the 
COMMON  blocks  will  be  given  in  Part  III. 

Subroutine  MATINV  is  a  general-purpose  matrix  inversion  subroutine 
borrowed  intact  from  the  computer  library  here  at  BRL.  Upon  return  from 
MATINV,  the  input  matrix  has  been  replaced  by  its  inverse. 

Subroutine  MERSO  is  a  general-purpose  numerical  integration 
subroutine  based  on  a  method  proposed  by  R.H.  Merson  of  Australia.  The 
method  is  a  fourth-order  member  of  the  Runge-Kutta  family,  requiring 
five  function  evaluations  at  each  integration  step.  The  subroutine 
adjusts  the  integration  step  size  automatically  to  obtain  a  predefined 
accuracy.  (All  of  this  is  transparent  to  the  FINLIE  user.) 

The  computer  library  at  BRL  contains  a  subroutine  MERSON  (see 
References  11  and  32  in  the  Bibliography)  for  performing  Merson 
integration.  Subroutine  MERSO  is  identical  to  MERSON  with  the  exception 
of  two  statements.  Firstly,  where  MERSON  has 

DIMENSION  T(IOO),  G(IOO),  S(IOO), 

MERSO  has  increased  the  three  dimensions  to  400  each.  Secondly,  where 
MERSON  has 


IF  (NT.LE.lOO)  GOTO  100, 

MERSO  compares  NT  with  400.  The  reason  for  these  changes  is  as  follows. 
The  size  of  the  three  arrays  T,  G  and  S  above  must  equal  or  exceed 

N5  5  N2  +  NA  (72) 


that  is,  the  number  of  differential  equations  (N2)  plus  the  number  of 
influence  equations  (NA) .  MERSON  requires  N5  <  100.  In  my  largest 
application  of  FINLIE  so  far,  N5  exceeded  100  (was,  in  fact,  368). 

Hence  the  minor  surgery  that  altered  MERSON  into  MERSO;  the  maximum 
permitted  value  of  N5  is  now  400.  This  value  of  400  appears  not  only 
in  MERSO  but  in  PARIS,  where  it  is  the  declared  dimensions  of  arrays 
U  and  DU  (see  the  Appendix).  Hence  the  user  can  change  the  upper  limit 
on  N5  by 


fi) 

and  (ii) 


changing  the  dimensions  of  U  and  DU  in  PARIS; 

changing  the  400  in  the  PARIS  statement  that  currently 
reads : 

IF(M5.LE.400)  GOTO  24 


and  (iii) 


changing  the  two  previously  mentioned  statements  in 
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MERSO: 


DIMENSION  T(400),  G(400),  S(400) 
and  IF(NT.LE.400)  GOTO  100 

Note  that  a  "nonstandard”  FORTRAN  function  appears  on  the  line  above 
statement  410  in  MERSO: 


H  =  SIGNl(H)  *  HMI 

Here  SIGNl  is  the  signum  function;  if  it  is  not  recognized  by  the 
user's  FORTRAN  compiler,  the  ahov v  statement  can  be  replaced  by 

IFCH.NE.O.)  H=5  GN(1.,H)*HMI 

The  use  of  multiple  arithmetic  and  logical  assignment  stai-ements  in 
MERSO  may  also  be  unacceptable  to  some  compilers.  In  a  multiple 
statement  of  the  form 


VN  =  =  V2  -  VI  =  expression, 


the  assignments  are  carried  out  from  right  to  left: 

VI  =  expression 
V2  *  VI,  etc. 

It  should  be  pointed  out  after  all  this  exposition  on  MERSO  that 
when  the  user  is  fitting  algebraic  or  transcendental  equations  rather 
than  differential  equations,  MERSO  is  not  needed  and  may  be  removed 
from  FINLIE. 


FINLIE  was  written  for  BRL's  CDC  computer  system  for  which  the 
single  piecision  of  real  numbers  is  approximately  14  decimal  digits. 

So  far,  this  has  proven  adequate  for  all  our  FINLIE  aoplications.  If 
the  user  is  working  with  a  machine  whose  single  precision  is  signifi¬ 
cantly  less  than  14  decimal  digits,  he  may  have  to  add  some  double¬ 
precision  declarations  to  the  subroutines  of  FINLIE.  One  source  of 
trouble  is  the  possibly  erratic  behavior  of  e  near  a  minimum,  due  mainly 
to  round-off  noise.  Hence  a  likely  candidate  for  double  precision 
is  array  GAMMA  in  subroutine  DUBLIN  (and  its  dummy  version  A  in 
subroutine  MATINV) .  A  more  complete  list  of  variables  that  may  require 
double  precision  includes: 


in  DUBLIN: 
in  LONDON; 
in  PARIS; 
in  MATINV; 


EA,  EB,  EPS,  GAMMA 
EP,  EPS,  RSQ 
RM,  RSQ 
A,  T1 


In  extreme  cases,  the  user  can  simply  double-precision  everything  in 
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bight;  this  may  be  inefficient  in  terms  of  storage,  but  it  could  save 
wear  and  tear  on  the  user. 


III.  OUTSIDE  FINLIE:  WHAT  THE  USER  MUST 
DO  FOR  FINLIE 

Assuming  that  the  FINLIE  user  is  given  a  set  of  equations  of  the 
form  (1)  or  (2)  -  or  equations  that  can  be  put  into  one  of  those  forms  - 
the  first  task  the  user  must  perform  is  to  derive  the  associated 
influence  equations,  as  indicated  in  parts  C  and  D  of  Section  II. 

The  next  task  is  to  write  a  FORTRAN  subroutine  defining  all  these 
equations  -  the  original  set  and  the  influence  equations  -  in  a  manner 
acceptable  to  FINLIE.  The  rules  for  constructing  this  subroutine  are 
slightly  different  for  sets  (1)  and  (2).  (We  assume  in  what  follows 
that  the  reader  has  some  familiarity  with  -  though  he  need  not  be  an 
expert  in  -  some  version  of  FORTRAN  equivalent  to  or  newer  than  FORTRAN 
IV.) 


A.  ROME:  The  User's  Subroutine  for  Fitting  Differential  Equations 

The  first  three  statements  of  ROME  have  the  form: 

SUBROUTINE  R0ME(N5,XE,U,DU) 

DIMENSION  U(N5),  DU (NS) 

COMMON/NAPLES/ PAR (40),  FLAG(60) 

It  should  be  noted  that  the  only  name  in  the  above  three  statements 
that  the  user  is  not  allowed  to  change  is  NAPLES.  All  other  names, 
including  ROME,  may  be  replaced  by  other  legal  FORTRAN  names  of  the 
user's  choice.  (Of  course,  the  distinction  between  integer  and  real 
names  should  be  maintained.) 

ROME  is  called  by  MERSO  (see  Figure  1)  and  hence  the  nature  of  the 
four  arguments  of  the  SUBROUTINE  ROME  statement  has  been  decreed  by 
MERSO.  The  first  three  arguments  are  input  to  ROME  (from  MERSO); 

N5  =  the  number  of  equations  to  be  defined  in  ROME;  N2 

(first-order  differential  equations)  plus  NA  (influence 
equations).  Thus  for  sample  set  (3),  the  value  of  N5 
is  2  +  8  =  10.  Note,  however,  that  this  argument  is  an 
integer  name,  not  an  integer  constant.  As  certain 
arrays  are  currently  dimensioned,  N5  cannot  exceed  400 
(see  the  pertinent  remarks  in  section  iI(M)). 


43 


'^^rKrmwniiiwi  iuii.jiiii».t»»  » 


XE  *  X,  the  independent  variable  value  at  which  the  N5 

equations  are  to  be  evaluated.  The  argument,  of  course, 
must  be  a  real  name,'  not  a  real  constant.  If  the 
independent  variable  does  not  appear  in  any  of  the 
equations,  then  argument  XE  will  not  be  used  in  the  body 
of  subroutine  ROME. 

U  ■  the  vector  of  N2  dependent  variables  and  NA  influence 
coefficients,  where 

Xj  =  U(J), 

Djk"  U(J  +  K*N2) 

[j=j-l,2,..N2  ^ 

K-k-l,2,..N23] 

Thus  for  sample  set  (3),  where  N2«>2  and  N23«4,  we  have 


U(1) 

= 

>^1 

U(2) 

= 

^2 

U(3) 

= 

"ll 

= 

»y,/9yio 

U(4) 

S 

^21 

= 

U(5) 

= 

^12 

= 

U(6) 

= 

D22 

8 

U(7) 

"l3 

8 

U(8) 

= 

*^23 

8 

3y2/3Ci 

U(9) 

= 

“l4 

X 

3yj/3c2 

U(10) 

s 

°24 

- 

3X2/902 

The  final  argument  of  ROME  is  an  output  (to  MERSO) : 

DU  =  the  derivative  vector  at  the  current  value  XE  of  the 
independent  variable,  where 

DU(J3  »  dU(J)/dx  (74) 

[J=1,2,..N5) 

Additional  input  to  ROME  comes  from  PARIS  via  the  labelled 
COMMON  block  NAPLES.  The  one  hundred  elements  of  the  NAPLES  block  are 
as  follows: 

PAR  =  a  vector  of  the  current  values  of  the  N3  parameters 
(not  paramics),  where  N3  <  4C.  For  sample  set  (3), 
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PAR(l)  = 

PAR (2)  = 

and  the  remaining  38  elements  of  PAR  are  undefined. 

FLAG  »  a  vector  of  N23  flags  (N23  <  60)  associated  with  the 
N23  single-round  paramic  set  Q,  Eq.  (12).  That  is, 
the  first  N2  elements  of  FLAG  are  associated  with  the 
N2  initial  conditions  and  the  remaining  N3  elements 
of  FLAG  are  associated  with  the  N3  parameters.  The 
value  of  FLAG(J)  is 

(i)  zero  if  the  value  of  the  corresponding 
paramic  is  fixed; 

or  (ii)  1.0  if  the  value  of  q.  is  to  be  adjusted 
by  the  fitting  proces^. 

For  sample  set  (3),  we  have 

FLAG(l)  =  flag  for 

FLAG (2)  =  flag  for 

FLAG (3)  =  flag  for  Cj 

FLAG (4)  =  flag  for  C2 

and  the  remaining  56  flags  are  undefined. 

Note  that  PAR  and  FLAG  are  inputs  to  ROME  from  FINLIE;  when 
writing  ROME,  the  user  assumes  that  the  two  arrays  already  contain 
their  proper  values.  In  the  case  of  the  initial  condition  flags, 
these  value:'  may  change  from  round  to  round.  For  example,  in  our 
tri-round  situation,  we  might  decide  to  make  a  computer  run  with  y^^^ 

for  round  El  and  y^^  for  round  E3  f ixed  at  specified  values.  Then 

FLAG(l)  0.0,  1.0,  1.0 
FLAG (2)  -  1.0,  1.0,  0.0 

for  rounds  El,  E2  and  E3,  respectively.  FINLIE  will  automatically 
change  the  values  of  FLAG(l)  and  FL.\G(2)  to  match  the  round  whose 
measured  data  is  currently  being  fitted.  Of  course,  FINLIE  can't 

fuess  what  the  user  wants  to  do;  it  must  be  told.  FINLIE  can  only 
efine  PAR  and  FLAG  on  the  basis  of  certain  inputs  given  to  it  by  the 
user's  main  program.  These  inputs  will  be  discussed  in  section  III(C). 

The  dimensions  of  PAR  and  FLAG  are  arbitrary  to  this  extent:  they 
can  be  changed  in  ROME  if  the  user  is  willing  to  make  all  the 
associated  changes  in  FINLIE.  To  save  space,  I  leave  the  nature  of 
such  changes  as  an  exercise  for  the  interested  reader.  The  simplest 


45 


course  is  to  make  no  changes  if  N2  <  40  and  N23  <  60. 

After  writing  down  the  first  three  statements  of  RCWE,  the  user 
is  ready  to  encode  the  body  of  the  subroutine:  the  statements  defining 
the  N5  elements  of  output  array  DU.  Consider,  for  example,  system 
(3).  For  convenience,  we  repeat  here  the  original  equations: 

Vl  “  1/72 

"  -Ci(l+C2y2)y2 

and  the  associated  influence  equations  (Table  III-A): 


(Dii)^ 

t: 

■‘^21/^2 

(D21)' 

= 

•c, (i+2c2y2)D2, 

CD, 2)' 

9 

■^22/^2 

(D22)' 

- 

-Ci(U2c2y23D22 

(D13)' 

=: 

■^23'^^'2 

(D,4)' 

s 

■^24/^2 

(D24)' 

* 

“C, (l+2c2y23D24  ~ 

For  these  equations,  a  likely  version  of  subroutine  ROME  is  given  in 
Table  VI. 
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Table  VI.  Subroutine  ROME  for  System  (3) 

SUBROm’INE  ROME  (N5,XE,U,DU) 

DIMENSION  l](N5),DUCN5) 

COMMON  /NAPLES/  Cl, C2, BLANK (38) , FLAG (60) 
V  =  U(2) 

DU(1)  =  l./V 
A1  =  C2*V 

DU(2)  =  -Cl*(l.  +  A1)*V 
A2  =  DU(1)**2 
A3  =  Cl*(l.  +  A1  +  Al) 

IF  (FLAG(l)  .EQ.  0.)  GOTO  10 
DU (3)  =  -A2*U(4) 

DU (4)  =  -A3*U(4) 

10  IF  (FLAG (2)  .EQ.  0.)  GOTO  20 
DU(5)  =  -A2*U(6) 

DU (6)  =  -A3*U(6) 

20  IF  (FLAG (3)  .EQ.  0.)  GOTO  30 
DU(7)  =  -A2*U(8) 

DU(8)  =  -A3*U(8)  -  (1.  +  A1)*V 
30  IF  (FLAG (4)  .EQ.  0.)  GOTO  40 
DU(9)  =  -A2*U(10) 

DU(10)=  -A3*U(10)  -  C1*V*V 
40  RETURN 
END 


Note  that  in  COMMON/NAPLES/  I  opted  to  write  the  forty- element 
parameter  set  in  the  form 


Cl,  C2,  BLANK(38) 

since  only  the  first  two  of  the  forty  elements  have  any  meaning.  I 
could  just  as  well  have  written  PAR(40)  in  the  COMMON  statement  and 
used  'r’AR(l)  and  PAR(2)  instead  of  Cl  and  C2  in  the  body  of  the 
subroutine.  Note  also  that  I  dimensioned  FLAG  as  60  even  though  the 
last  56  elements  are  meaningless.  This  was  a  courtesy  to  our  CDC 
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FORTRAN  compiler,  which  likes  all  COMMON  blocks  of  the  same  name 
(NAPLES  in  this  case)  to  have  the  same  length.  The  compiler  doesn't 
insist  when  you  break  this  rule,  but  it  comments  on  your  bad  form. 

Because  ROME  will  be  called  many  times  by  MERSO  during  the  course 
of  the  numerical  integration,  the  user  should  take  the  time  to  make 
ROME  as  efficient  as  practicable.  For  large  and  labyrinthian  systems 
of  equations,  a  worthy  ROME  isn't  built  in  a  day. 

One  of  the  aids  to  efficiency  in  ROME  is  the  FLAG  vector.  Note 
that  if  any  paramic  value  is  fixed  during  a  computer  run  (that  is,  if 
the  associated  flag  value  is  zero),  the  influence  equations  for  that 
paramic  need  not  be  calculated.  Hence  the  FLAG  vector  can  -  and  in  my 
opinion  should  -  be  used  as  indicated  in  Table  VI  to  avoid  these 
unnecessary  calculations.  The  general  rule  is  that  if  FLAG(J)  is  zero, 
then  ROME  need  not  evaluate  DU (LA)  through  DU (LB),  where 

LA  =  (J  X  N2)  +  1 

LB  =  (J  X  N2)  +  N2  =  LA  +  (N2-1) . 

Of  course,  if  the  user  is  convinced  that  he  will  never,  ever  want  to 
hold  fixed  the  value  of  some  paramic,  he  can  omit  the  corresponding 
IF- statement  from  ROME. 

Some  systems  of  equations  may  involve  constants  whose  values  are 
always  fixed  (that  is,  never  adjusted  by  FINLIE)  and  yet  these  values 
may  change  from  run  to  run.  It  would  be  possible  -  but  not  too  bright  - 
to  handle  sucb  i  constant  as  a  fixed  parameter:  a  parameter  whose 
associated  flag  is  always  zero.  A  better  approach  is  to  pass  any  such 
constant  directly  from  the  user's  main  program  to  ROME  through  a  new 
labelled  COMMON  block  (see  block  CAIRO  in  Figure  1).  Of  course,  if 
a  constant  will  never  change  from  run  to  run,  it  need  only  be  defined 
within  ROME. 

A  final,  rather  minor  comment:  Sample  set  (3)  is  one  of  those 
cases  where  the  input  argument  XE  is  not  used  in  the  body  of  subroutine 
ROME,  simply  because  the  independent  variable  does  not  appear  explicitly 
in  the  N5  equations  of  this  example. 


B.  ROMA:  the  User's  Subroutine  for  Fitting  Algebraic  or  Transcendental 
Equations 

Many  of  the  comments  in  the  previous  section  concerning  ROME  apply 
to  ROMA  as  well.  Hence,  if  the  reader  has  skipped  over  that  section 
because  his  interest  in  fitting  differential  equations  is  minimal, 
he  may  have  missed  something  noteworthy.  Or  possibly  not. 

The  first  three  statements  of  ROMA  have  the  form: 
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SUBROUTINE  ROMA  (COND,XO, XE,U) 

DIMENSION  COND(n2),  U(n5) 

COMMON/NAPLES/PAR (40),  FLAG(60) 

The  first  three  arguments  in  the  SUBROUTINE  statement  are  inputs 
(from  PARIS) : 

COND  =  a  vector  of  N2  current  initial  condition  values.  For 
sample  set  (4), 

COND(l)  = 

COND (2)  = 

For  multi-round  situations,  the  initial  conditions 
change  with  the  round  as  well  as  with  the  current  state 
of  the  fitting  process.  FINLIE  supplies  the  proper 
COND  vector  to  ROMA  automatically. 

XO  =  Xq,  the  independent  variable  value  at  which  the  initial 
conditions  apply.  (This  one  value  must  apply  to  all 
rounds. ) 

XE  =  X,  the  independent  variable  /alue  at  which  the  equations 
are  to  be  evaluated. 

The  final  argument,  U,  is  an  output  vector  defined  exactly  as  in  the 
previous  section  for  subroutine  R(X4E. 

In  the  DIMENSION  statement,  n2  and  nS  denote  the  values  of  N2  and  N5, 
respectively.  (Actually,  on  the  CDC  system  and  on  most  other  computers, 
a  one-dimensional  argument  array  in  a  subroutine  need  not  be  declared 
at  its  maximum  size;  the  value  1  is  adequate.) 

The  labelled  COMMON  block  NAPLES  brings  to  ROMA  the  arrays  PAR 
and  FLAG,  defined  in  the  previous  section. 

The  body  of  subroutine  ROMA  consists  of  the  statements  defining 
the  needed  elements  of  array  U.  Consider,  for  example,  system  (4). 

For  convenience  we  repeat  here  the  original  equations  (4)  and  the 
needed  influence  equations  (Table  111(B)): 

XlO  -  (b/Cj)(u-l) 

(bu  -  c^)"^ 

1. 


11 


12 


2 

Dj3  =  (b/Cj  ) [1-u+CjCx-Xq)u] 

®14  “  (u-l)/Cj,  -  (x-Xq) 

where 

u  =  exp  [Cj(x-Xq)] 

b  =  Cy^o)"*  »  <=2 

For  these  equations,  a  likely  version  of  subroutine  ROMA  is  given  in 
Table  VII. 


Table  VII.  Subroutine  ROMA  for  System  (4) 


SUBROUTINE  ROMA  (COND, XO, XE,U) 
DIMENSION  COND(2),U(10) 


COMMON  /NAPLES/  C1,C2, BLANK(38) , FLAGC60) 


TO 

=  COND(l) 

VO 

=  COND (2) 

XA 

=  XE  -  XO 

z 

=  EXP(C1* 

XA) 

B 

=  C2  +  1. 

/VO 

U(l) 

U(2) 

IF 

(FLAGd) 

.NE.  0.) 

U(3) 

IF 

(FLAG (2) 

.NE.  0.) 

U(5) 

IF 

(FLAG (3) 

.NE.  0.) 

U(7) 

IF 

(FLAG (4) 

.NE.  0.) 

U(9) 

RETURN 

END 

TO  -  C2*XA  +  B*(Z  -  l.)/Cl 
l./(B*Z  -  C2) 

1. 

(1.  -  Z)/(C1*V0**2J 

B*(l.  -  Z  +  C1*XA*Z)/(C1**2) 

(Z  -  l.)/Cl  -  XA 


As  discu.ssed  in  section  11(D),  FINLIE  does  not  require  expressions 

for  the  influence  coefficients  D.,  when  j  is  greater  than  Nl.  Hence 

J  K 

in  this  sample  ROMA,  where  Nl=l,  the  D2j^  equations  (namely,  the  equations 
for  U(4),  U(6),  U(8)  and  U(10))  are  simply  omitted  from  the  subroutine. 


As  with  ROME  in  the  previous  section,  the  FLAG  array  in  ROMA  is 
used  to  avoid  calculating  D  when  the  value  of  paramic  q,  is  fixed, 

J  K  K 

Also  as  with  ROME,  any  needed  "changeable  constants"  can  be  passed 
directly  from  the  user's  main  program  to  ROMA  through,  say,  the 


labelled  COMMON  block  CAIRO. 
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C.  Calling  Subroutine  DUBLIN 

kfteT  the  user  has  written  his  subroutine  defining  the  equations 
to  be  fitted,  his  next  task  i  s  to  write  a  program  unit  •-  we  assime  a 
main  program  -  that  utilizes  FINLIE.  Before  discussing  this  main 
program  as  a  whole,  we  will  concentrate  on  one  statement  within  that 
main  program:  the  CALL  DUBLIN  statement. 

This  statement  is  the  link  between  the  user  and  FINLIE.  It  can  be 
written  in  the  form 

CALL  DUBLIN  (R0ME,NF,N1,N2,N3,N7,N8,NR,NM,X0, 
X,Y,F,NW,W,P,RL,NC,YC,R,RS,EPS, 

SIG,EK,NS) 

where  all  integer  names  happen  to  start  with  the  letter  N.  The  first 
fourteen  of  the  twenty-five  arguments  are  inputs. 

[1]  ROME  is  the  name  of  the  subroutine  [written  by  the  user)  that 

defines  the  equations  to  be  fitted  (See  sections 
IIT(A-B)). 

The  values  of  the  remaining  thirteen  input  argiments  must  be  established 
in  the  user's  main  program  before  DUPLIN  is  called.  These  values  will 
not  be  changed  by  FINLIE;  hence  actual  values  rather  than  names  may  be 
used  for  arguments  [2]  through  [8],  [10]  and  [14]  below. 

[2]  ^  is  a  flag  that  indicates  the  nature  of  the  equations 

to  be  fitted: 

NF=0  if  the  fitting  equations  are  algebraic  or 
transcendental  (System  (2)); 

NF=1  if  the  fitting  equations  are  differential 
equations  (System  (1)). 

[3]  ^  is  the  number  of  measured  dependent  variables  in  the 

system,  where 

1  <  N1  <  10  (75) 

(The  upper  bound  on  N1  -  and  the  upper  bounds  indicated 
for  some  of  the  other  arguments  defined  below  -  can  be 
increased  only  by  delving  into  FINLIE.)  N1  must  have 
the  same  value  for  each  round;  FINLIE  insists  that  the 
same  dependent  variables  be  measured  for  each  round  used 
in  the  fitting  process. 

[4]  N2  is  the  total  number  of  dependent  variables  in  the  system, 

where 
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N1  <  N2  <  20  (76) 

[5]  ^  is  the  naxlmum  number  of  parameters  (not  paramics)  whose 

values  can  be  determined  from  the  fit,  where 

0  <  N3  <  40  (77) 

I  use  the  word  "maximum"  above  because  the  actual  number 
of  parameters  to  be  determined  in  the  course  of  a  computer 
run  may  be  less  than  N3.  The  user  specifies  (by  argument 
F,  to  be  discussed  below)  which,  if  any,  of  the  parameters 
and  initial  conditions  are  to  be  held  fixed  at  their  input 
values  and  which  are  to  be  adjusted  by  FINLIE  during  the 
run.  Input  N3  is  the  total  number  of  parameters:  those  to 
be  adjusted  plus  those  held  fixed.  (If  input  N3  is  zero  - 
the  lower  limit  in  inequality  (77)  -  then  presumably 
there  is  at  least  one  initial  condition  to  be  determined; 
otherwise  there  would  be  r.o  reason  for  running  the 
program. ) 

[6]  ^  is  the  number  of  roivs  declared  in  the  user's  main 

program  for  the  two-dimensional  arrays  Y,  W  and  R  defined 
below  as  arguments  [12],  [15]  and  [20],  respectively. 

As  we  will  uce,  these  three  arrays  serve  as  Nl  by  N4 
matrices.  At  first  glance,  then,  it  might  seem  that 
N7=N1.  However,  the  user  may  not  want  to  restrict  his 
main  program  DIMENSION  statement  to  the  current  values 
of  Nl  and  N4.  It  is  often  more  convenient  to  dimension 
arrays  at  their  largest  anticipated  sizes.  For  example, 
in  our  recurring  case  where  Nl=l  and  N4=16,  the  user 
might  want  to  dimension  arrays  Y,  K  and  R  as,  say  (2,50) 
rather  than  (1,16).  FINLIE  will  go  along  with  this  sort 
of  thing,  but  it  wants  to  be  told  about  it.  Thus  if  the 
user  dimensions  Y,  W  and  R  as  (2,50),  he  must  set  N7 
eaual  to  2.  In  general,  then,  N7  >  Nl .  (The  declared 
column  dimension  for  the  three  arrays  -  say,  50  -  is  not 
needed  by  FINLIE.  The  declared  row  dimension  is 
sufficient  -  assuming  the  computer  stores  matrices  in  the 
usual  way,  that  is,  by  columns  -  to  maintain  notational 
row-column  agreement  between  calling  program  and 
subroutine.  Neither  is  FINLIE  interested  in  the  declared 
dimensions  of  its  vector  arguments.*) 


*fn  FTNLIEj  I  have  declared  1  as  the  last  (right-most )  di-mension  of 
subroutine  durtny  arg’ument  arrays.  This  is  fairly  oormon  FORTRAN  4 
practice,  but  FORTRAN  77  prefers  an  asteid.sk:  Y(N7,  a)  instead  of 
Y(N?,1). 
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[7]  1^  is  the  number  of  rows  declared  in  the  user's  main 

program  for  the  two-dimensional  array  YC  defined 
below  as  argument  [19].  Array  YC  serves  as  an 
N2  by  N4  matrix;  hence  N8  >  N2.  (See  the  comments 
for  argument  [6]  above.) 

[8]  is  the  number  of  data  rounds  to  be  considered 
simultaneously,  where 

1  <  NR  <  (60-N?;)/N2  (78] 

The  right-half  of  this  double  inequality  may  seem  a 
rather  strange  condition  to  spring  upon  the  reader. 

Until  now,  no  limit  has  been  implied  on  the  number  of 
rounds.  The  basic  condition  (somewhat  concealed  in 
(78))  is 

N  <  60  (79) 

where  N  is  the  total  number  of  paramics,  (NR  x  N2)  +  N3. 
Condition  (79),  like  the  limits  on  Nl,  N2  and  N3,  is  a 
result  of  arbitrary  DIMENSION  decisions  that  had  to  be 
made  when  constructing  FINLIE.  Since  N  is  not  itself 
an  input  to  DUBLIN,  I  have  simply  converted  (79)  to  the 
equivalent  form  (78).  By  satisfying  (78),  the  user  can 
be  sure  that  (79)  is  also  satisfied.  For  saxriple  set  (3) 
or  (4)  we  must  have  NR  <  (60-2) /2=29  For  the 

associated  data  of  Table  I,  we  have  NR=3,  well  below  the 
maximum  permitted.  Recall  that  the  data  for  an 
individual  round  solely  determine  the  initial  conditions 
for  that  round,  but  combine  with  the  data  from  all  the  other 
rounds  to  determine  the  parameters. 

[9]  ^  is  a  vector  of  NR  elements,  where 

NM(J)  =  the  number  of  data  points  (that  is,  the 
number  of  independent  variable  values  at 

which  measurements  were  taken  )  foi  the  J-th 
round . 

Thus  for  the  sample  data  of  Table  I,  the  user's  main 
program  must  set 

NM(1)  =  5 
NM(2)  =  6 
NM(3)  =  5 

FINLIE  determines  N4,  the  total  number  of  data  points, 
by  summing  the  NM  components: 
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N4 


(80) 


NM(J) 


The  user  must  insure  that  N4  satisfies  the  inequalities 


N  <  N4  <  1000 
N4xN5  <  10000 


J 


and  that 

N1  X  (MAX. ELEMENT  OF  NM)  <200  )  , 

N2  X  (MAX.  ELEMENT  OF  MM)  <  400  /  '• 

Again,  these  restrictions  are  the  result  of  arbitrary 
DIMENSION  statements  in  FINLIE. 


XO  is  the  independent  variable  point  x«  at  which  all 

initial  conditions  apply.  The  same^x^  must  apply  to 
all  rounds. 


)(  is  a  vector  of  the  N4  independent  variable  values 
x^  at  which  measurements  were  taken.  The  first  NM(1) 

values  in  X  are  the  first-round  values,  in  increasing 
order; 


X(M-l)  <  X(M),  M  =  2,3,...NM(1) 

The  next  NM(2)  values  of  X  are  the  second-round  values, 
also  in  increasing  order  among  themselves: 

X(M-l)  <  X(M),  M  -  NM(1)  =  2,3,...NM(2) 

and  so  on.  For  the  Table  I  data,  we  have 

X(M)  X,  m=l,2,...16. 
m  ’ 

is  an  N1  by  N4  matrix  of  measured  dependent  variable 
values,  where 

Y(T,M)  =  the  measured  value  of  y^  at  X(M) 

For  the  Table  I  data,  we  have  Y(1,M)  =  m=l,2,...16. 

£  is  a  vector  of  N  flags  associated  with  the  paramic  point 
P  (argument  [16]  below),  where 


F(J)  =  0.0  if  the  input  value  of  P(J)  is  to 
be  held  fixed; 

=  1.0  if  the  input  value  of  P(J)  is  to  be 
adjusted  by  the  fitting  process. 

[14]  M  is  £  weight  flag  associated  with  matrix  W  (argument  [15] 
below) . 


NW  =  0  if  the  user*?  weights  (already  stored 
in  matrix  W)  are  to  be  used  by  FINLIE; 

=  1  if  all  weights  are  unity  (in  which  case, 
the  user  need  not  store  1.0's  in 
matrix  W  before  calling  DUBLIN). 

Recall  the  comments  in  the  vicinity  of  Eqs.  (8)  and  (9') 
regarding  weights.  The  important  point  is  that  the 
"easy"  way  out  -  assigning  unit  weights,  merely  by 
setting  NW=1  -  will  often  lead  to  a  poor  fit.  Give 
some  minimum  consideration  to  the  possibility  of 
unequal  uncertainties  in  the  measurements,  particularly 
when  more  than  one  variable  has  been  measured  (N1  >  1) . 

Tlie  fifteenth  argument  of  DUBLIN  may  or  may  not  be  defined  by  the 
user  before  DUBLIN  is  called: 

[15]  W  is  the  N1  by  N4  matrix  of  weights  associated  with  input 

Y  (argument  [12]  above).  The  user  has  a  choice  to 
make.  If  each  of  the  N1  by  N4  measurements  in  matrix 

Y  can  be  assigned  unit  weight;  that  is,  if 

W(I,J)  =  1.0 

then  the  user  need  not  define  the  W  array.  Simply  set 
NW  (argument  [14]  above)  to  1.  If,  on  the  other  hand, 
the  user  decides  that  one  or  more  of  the  weights  must 
differ  fi'om  1,  then  the  user  must  define  the  entire 
array,  subject  to  the  conditions  that  each  weight  be 
nonnegative  and  that 

[W(I,J)]j  -  [l./Y(l,.J)^]j 

See  the  comments  near  Eqs.  (8)  and  (9). 

The  next  three  arguments  of  DUBLIN  are  input/output.  That  is,  the 
user  must  define  them  before  the  first  CALL  DUBLIN  statement,  but 
FINLIE  will  change  their  values. 
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[16]  P  is  the  current  N-dimensional  paramic  point  P,  Eq.  (5), 


where 

P(l],  .  p(N2)  =  first  round  initial 

conditions, 

P(l+N2), . P(2*N2)  =  second  round  initial 

conditions, 


P(1+(NR-I)*N2)  .  .  . 

P(NR*N2)  = 

last  round  initial 
conditions. 

P(1+NR*N2)  . 

P(N) 

the  parameters. 

Thus  for  sample  system  (3)  or  (4)  and  for  tri-round  data, 
the  eight  elements  of  P  are  given  by  Eq.  (10).  Clearly, 

P  is  the  essential  argument  in  the  CALL  statement;  the 
other  argioments  play  a  necessary  hut  supportive  role. 

As  indicated  in  Part  I  of  this  report,  an  effort  should  be 
made  to  find  suitable  starting  values  for  the  elements 
of  P.  Not  all  first  estimates  will  lead  to  the  right 
answer.  Each  time  the  program  returns  fiom  DUBLIN, 
array  P  will  contain  an  updated  point.  More  precisely, 
the  first  call  to  DUBLIN  is  a  special  situation  and  P  is 
unchanged  upon  return.  Thereafter,  each  call  serves  to 
update  P.  For  more  on  this  first  call,  see  argument 
[18]  below.  In  general,  then,  each  DUBLIN  call  after  the 
first  advances  P  one  step  on  the  road  to  the  solution. 
■DOTCTK  should  be  called  repeatedly  (say,  in  a  DO-loop) 
until  convergence  is  achieved.  Not  all  elements  of  P 
will  necessarily  change  with  the  iteration.  If  input 
F(J)  is  zero  (see  argument  [13]  above),  then  the  original, 
user-assigned  value  of  P(J)  will  be  maintained  no  matter 
how  many  times  DUBLIN  is  called, 

[17]  ^  is  a  Marquardt  argument.  Before  DUBLIN  is  called  the 

first  time, 

(i)  Set  RL  -  0.0  if  the  Marquardt  algorithm  is  to  be 
omitted  from  the  fitting  process 
(that  is,  if  the  user  wants  FINLIE 
to  fit  by  differential  corrections, 

Eq.  (56),  rather  than  by  Marquardt 
interpolacion,  Eq.  (58)).  In  this 
case,  RL  will  remain  at  zero. 

(ii)  Set  RL  •-=  1.0  if  the  Marquardt  algorithm  is  to  be 

used.  Upon  the  first  return  from 
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DUBLIN,  RL  will  have  the  "starting" 

X  value  of  0,01.  (On  subsequent 
calls  to  DUBLIN,  the  input  val'  e  of 
RL  is  inunediately  divided  by  ten; 
hence  the  true  starting  vs.ue  of  X 
is  0.001,  as  indicated  in  Eq.  (61)). 
Upon  the  second  and  subsequent 
returns  from  DUBLIN,  RL  will  have  the 
value  of  A  used  to  obtain  the  point 
returned  in  array  P.  Note  that 
since  FINLIL  changes  RL,  a  name  (not 
the  value  1.0)  must  be  used  in”  the 
CALL  list. 

[18]  is  a  "first  call"  flag.  The  user  must  set  NC=0  initially. 
This  value  alerts  FINLIE  to  the  fact  that  it  is  being 
called  for  the  first  time.  FINLIE  behaves  differently 

on  this  first  call  than  it  does  on  all  subsequent  calls. 

In  particular,  on  the  first,  call,  FINLIE 

(i)  sets  all  elements  of  argument  W  to  unity  if 
NW=1; 

(ii)  sets  argument  RL  to  0.01  if  the  input  RL  is  1.0 

(that  is,  if  Marquardt's  algorithm  i.s  to  be  used); 

(iii)  determines  the  number  of  paramics  to  be  adjusted 
(the  total  number  minus  the  number  of  parajnics 
held  fixed)  and  stores  this  value  back  in 
argument  NC  (hence  use  a  name,  not  the  integer 
zero,  for  the  "first  call"  flag  in  the  CALL  list); 

(iv)  evaluates  the  next  five  argiiments  in  the  CALL 
list  (YC,R,RP,EPS  and  SIG,  all  described  below) 
at  the  input  point  P^. 

Note  that  FINLIE  does  not  update  the  input  point  P^^  on 

this  first  call:  P^  goes  in  and  comes  back.  The 

paramics  are  updated  only  on  the  second  and  subsequent 
calls. 

The  remaining  seven  arguments  of  DUBLIN  are  outputs,  evaluated  at 
the  current  value  of  P. 

[19]  ^  is  an  NP  by  N4  matrix  of  computed  dependent  variable 

values,  based  on  the  current  point  P,  where 

YC(J,M)  =  the  computed  value  of  y^  at  X(M) 

Thus  for  the  Table  I  data  in  our  examples. 
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Yca,MD  -  Xj  (X^,P) 

YC(2,M]  -  y2 

[m  =  1,2, . . .16] 

When  fitting  differential  equations,  FINIIE  obtains  the 
YC  values  by  numerical  integration  of  system  (1] ;  when 
fitting  algebraic  or  transcendental  eqi.ations,  FINLIE 
obtains  YC  directly  from  the  equation  set  (2). 

[20]  ^  is  an  N1  by  N4  matrix  of  residuals  of  the  fit,  where 

RCI.M)  =  YCI.M)  -  YC(I,M]  (831 


[21]  ^  is  a  vector  of  N1  nondimensional  error  measures 

associated  with  the  N1  measured  dependent  variables, 
where 


RS(I] 


=  that  parr  of  e  (see  Eq.  (7)  and  argument 

[22]  below)  that  can  be  attributed  to  the 
fit  on  y^ 


N4 


E 

M=1 


W(I,M)  [Ra,M)l^ 


(.84) 


[22]  EPS 


is  e(P],  the  nondimensional  sum  of  the 
of  the  residuals  of  the  fit  (Eq.  (711, 


weighted  squares 
where 


EPS 


N1 

Rs(i) 


(8S) 


If  the  Marquardt  feature  is  being  used  (see  argument 
[17]),  then  after  the  first  call,  DUBLIN  should  return  an 
EPS  no  greater  than  the  input  EPS. 


[23]  SIG 


[24]  EK 


is  the  estiir>ated  standard  deviation  of  the  fit 
(Eq.  (69)),  where 


SIG 


EPS 


] 


I 


'2 


(86) 


is  a  vector  of  crude  estimates  of  the  errors  in  the  N 
paramics  of  point  P,  where 


EK(K)  =  the  estimated  standard  deviation 
in  paramic  P(K) 
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as  defined  in  Eq.  (71) 


[2ri]  is  a  convergence  flag.  Before  returning  to  the  user's 
main  program,  FINLIE  will  set 

I1S=0  if  the  process  has  not  yet  converged  by  criterion 
(68),  but  there  is  still  hope.  FINLIE  is  saying 
in  effect,  "Nothing  obvious  has  gone  wrong  yet 
so  give  DUBLIN  another  call." 

NS-1  if  all  output  arguments  (except  this  one)  are 
invalid.  Usually  this  liappens  when  some  input 
argument  is  invalid.  (FINLIE  performs  a  few 
simple  checks  to  spot  invalid  inputs.)  If 
DUBLIN  returns  an  NS  value  of  1,  the  main  program 
should  take  some  special  action  (e.g.,  STOP). 

NS=2  if  the  latest  iteration  has  satisfied  convergence 
criterion  (68).  If  the  user  is  willing  to  accept 
this  criterion,  his  main  program  should  step 
calling  DUBLIN  when  NS=2.  If  the  user  is 
imposing  some  more  stringent  convergence 
criterion  of  his  own,  he  should  regard  NS=2  as 
having  the  same  meaning  as  NS=0. 

To  illustrate  the  use  of  these  twenty-five  arguments,  consider  our 
sample  systems  (3)  and  (4)  with  three-round  data.  Assume  that  in  the 
calling  program,  arrays  Y,  W,  R,  and  YC  have  been  dimensioned  as 
(2,50).  Then  for  system  (3)  and  subroutine  ROME,  we  can  write 

CALL  DUBLIN(ROi'iE,l,l,2,2,2,2,3,NM,0.0,X,Y,F, )  ,W,P.RL,NC, 
YC,R,RS,EPS,SIG,EK,NS) 

For  system  (4)  and  subroutine  RCWA,  only  the  firjt  two  arguments  above 
are  changed: 


CALL  DUBLIN (ROMA, 0,. . .) 


D,  Writing  the  Program  that  Calls  DUBLIN 

In  this  final  section,  a  typical  main  program  for  utilizing  FINLIE 
is  broken  down  into  six  steps.  Some  of  these  steps  are  essential, 
others  are  optional. 

Step  (1).  Dimension  all  ten  arrays  appearing  3n  the  CALL  DUBLIN 
statement: 


Dir-iENSTON  NM(nr),  X(n4),  Y(nl,n4),  F(n),  W(r.l,n4), 
P(n),  YC(n2,n4j,  R(nl,n4),  RS(nl),  EK(n) 
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where  small-letter  dimensions  above  denote  constants  no  less  than  the 
values  of  the  corresponding  capital- letter  names.  That  is,  nr  >  NR, 
etc.  As  we  have  mentioned,  it  is  often  useful  to  dimension  arrays  larger 
than  their  current  working  sizes.  For  example,  in  our  tri-round  test 
cases  (S)  and  (4),  we  might  write: 

DIMENSION  NMCS),  X(S0),  Y(2,50),  F(8),  W(2,50),  P(8),  YC(2,50), 
RC2.50),  RS(2),  EK(S) 

This  would  allow  for  up  to  five  rounds  (nr«5),  f if ty  measurement  points 
(n4»50)  and  two  measured  variables  (nl"2).  Note  that  the  values  given  to 
the  row  sizes  nl  and  n2  in  this  DIMBNSI014  statement  become  the  values 
of  arguments  N7  and  N8  when  DUBLIN  in  called.  On  the  other  hand,  the 
dimensions  allotted  above  to  the  vector  arguments  and  to  the  columns 
of  the  matrix  arguments  are  of  no  interest  to  FINLIE. 

Step  (2) .  Declare  in  an  EXTERNAL  statement  the  user  subroutine 
whose  name  will  be  passed  to  DUBLIN.  Thus  for  sample  set  (3)  and  the 
corresponding  ROME  (Table  VI) ,  we  would  write 

EXTERNAL  ROME 

and  similarly  for  set  (4)  and  RCWA. 

Step  (3) .  Establish  initial  values  for  seven  DUBLIN  arguments: 

NM,  X,Y,F,P,RL,NC 

and  if  necessary,  for  an  eighth  argument:  W.  There  is  no  standard 
coding  for  obtaining  the  values  of  these  arguments;  the  technique  will 
vary  with  the  situation.  For  example,  initial  estimates  for  array  P 
might  be  read  in  at  this  .stage,  or  they  might  be  obtained  by  calling 
some  subroutine  whose  sole  purpose  is  to  derive  adequate  estimates 
from  the  data.  For  simplicity,  let's  assume  that  in  our  main  program 
for  sample  set  (3)  or  (4), 

(a)  the  arrays  NM,X,Y,F  and  P  are  read  in; 

lb)  RL  and  NC  are  defined  explicitly: 

RL  =  1.0 
NC  =  0 

(c)  array  W  is  not  defined  (since  argument  NW  will  be  1  in 
the  CALL  statement) . 

Note  that  the  values  of  the  remaining  nine  input  arguments: 

NF,  Nl,  N2,  N3,  N7,  N8,  NR,  XO,  NW 
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can  be  established  in  the  CALL  DUBLIN  statement  itself. 


Step  (4) ,  Write  column  headings  for  everything  of  interest  that  will 
be  deteralned  at  the  end  of  each  iteration.  Of  course,  "interest"  is 
subjective.  One  user  may  want  a  detailed  print-out  of  the  progress 

A 

from  Pq  to  P;  a  less  inquisitive  user  may  care  only  for  what  pertains  to 

the  final,  converged  point.  Personally,  for  each  iteration,  I  like  to 
see: 


(a)  the  iteration  niomber  i  (i=0, 1, 2, . . . ) 

(b)  the  N  elements  of  point  P^ 

(c)  the  value  of  Marquardt’s  X  required  to  produce 

(d)  the  residual  function  e(P.)  and/or  the  standard  deviation 

of  the  fit  s(Pj^).  ^ 

These  desiderata,  then,  determine  my  column  headings.  (Of  course, 
if  what  I  want  to  see  can  not  be  conveniently  spread  across  a  single 
output  page,  then  some  of  the  results  of  each  iteration  have  to  be 
saved  -  by  storing  them  in  additional  arrays  -  so  that  they  can  be  printed 
later  on  a  second  page.) 

Step  (5) .  Program  the  DO- loop  that  calls  DUBLIN.  For  our  sample 
set  (3) ,  we  might  write: 

DO  60  K=l,26 

CALL  DUBLIN  (ROME, 1 , 1 , 2, 2, 2, 2, 3,NM, 0. 0, X, Y,F, 

1  1,W,P,RL,NC,YC,R,RS,EPS,SIG,EK,NS) 

NP0INT=K-1 

WRITE (6, 100) NPOINT , P , RL , SIG 
IF  (NS-1)  60,70,80 
60  CONTINUE 


WRITE  (6,101) 

C - The  above  is  a  warning  that  the  process  has 

C - failed  to  converge  in  25  iterations. 

GOTO  80 

70  WRITE  (6,102) 

C  - The  above  is  a  warning  that  something  is  wrong. 

STOP 


80  CONTINUE 

In  the  above  code,  DUBLIN  will  be  called  until  output  argument  NS 
equals  1  or  2,  or  until  the  DO-loop  variable  K  exceeds  26,  whichever 
occurs  first.  (The  limit  26  -  that  is,  25  iterations  -  is  arbitrary; 

1  is  not  enough,  lO^^  is  too  many.)  After  each  iteration,  we  obtain  a 
print-out  -  presumably  under  the  proper  column  headings  -  of  NPOINT 
(the  number  of  iterations),  the  N  elements  of  the  current  point  P,  and 
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finally  the  X  and  s  values  at  the  current  point.  The  first  line  of 
this  print-out,  where  NPOINT«0,  gives  the  initial  estimated  values  of 
the  paramics.  If  NvSal  at  the  end  of  any  iteration,  the  program  stops; 
otherwise  the  program  moves  eventually  to  statement  80.  Note;  the 
principal  results  of  the  fitting  process  are  the  final  printed  values 
of  the  N  paramics.  All  else  is  in  a  sense  window-dressing. 

Step  (6) .  Write  anything  else  of  Interest.  My  usual  scheme  is  as 
follows ; 

(a)  aligned  under  the  final  paramic  values  (but  with  one 
line  skipped  for  clarity),  I  write  the  corresiponding  crude  error 
estimates  contained  in  array  EK.  (If  flag  array  F  is  of  interest,  the 
elements  of  F  can  be  written  on  the  next  line,  again  aligned  under  the 
corresponding  paramic  values.) 

(b)  on  a  new  output  page,  I  write  the  values  stored  in  arrays 

X,  Y,  YC,  R  and  (possibly)  W, 

one  line  for  each  X  value.  (In  multi-round  fits,  I  skip  a  line  for 
clarity  at  the  end  of  the  data  for  each  round.) 

(c)  wherever  convenient,  I  write  the  suitably  labelled  values 
of  some  or  all  of  the  following; 

N1,N2,N3,N4,N,NC,NM,NR,NW,RS,X0 
IV.  SUMMARY 

The  recent  patter  of  tiny  details  has  very  likely  blurred  the  big 
picture.  To  review,  then,  assume  that  the  reader  has  a  problem 
reducible  to  fitting  a  set  of  equations  of  the  form  (1)  or  (2)  to 
measured  data.  Further  assume  that  this  reader--an  adventurous  spirit- - 
decides  to  use  FINLIE  to  solve  the  problem.  Then  this  invoker  of  FINLIE 
must ; 

(a)  derive  the  related  set  of  influence  equations  (Section  II,  C 
or  D); 

(b)  write  a  FORTRAN  subroutine  that  lists  the  original  equations 
and  the  related  influence  equations  (Section  III,  A  or  B) ; 

(c)  write  a  FORTRAN  main  program  (Section  III,  D)  that  will; 

(i)  furnish  adequate  initial  estimates  of  the  parameters  and 
initial  conditions; 

(ii)  specify  which,  if  any,  of  these  estimates  are  to  be 
adjusted  by  FINLIE; 

(iii)  assign  weights  to  the  measurements  (if  the  weights  are 
not  all  equal) ; 

(iv)  call  subroutine  DUBLIN  (Section  III,  C)  in  a  DO-loop; 
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(d)  submit  the  entire  program  (main,  FINLIE  and  equation-defining 
subroutine)  to  the  computer  and  ponder  the  ensuing  output. 

This  output  will  take  one  of  four  forms,  listed  in  decreasing 
order  of  desirability; 

(1)  convergence  to  the  right  answer; 

(2)  failure  to  converge  in  the  specified  number  of  itera¬ 
tions  (sometimes  achieving  an  apparent  oscillation  about  an  answer); 

(3)  divergence  (the  program  crashes) ; 

(4)  convergence  to  the  wrong  ansv:er. 

Result  (1)  above- -convergence  to  the  right  answer- -should  prevail 
when  all  of  the  following  hold; 

(a)  the  measured  data  are  a  representative  sample  of  the 
total  behavior  they  are  meant  to  define.  (An  elementary  violation 
would  be  measurements  taken  every  t  seconds  on  a  periodic  variable  of 
period  t.) 

(b)  the  measured  data  are  free  of  gross  errors. 

(c)  "least  squares"  is  a  suitable  fitting  criterion.  (This 
implies  that  the  measurement  errors  possess  certain  statistical  traits; 
however,  the  degree  to  which  the  errors  must  possess  these  traits  in 
order  to  be  considered  amenable  to  least  squares  is  a  matter  of  judg¬ 
ment  . ) 

(d)  the  fitting  equations  with  their  associated  parameters 
are  appropriate  for  describing  the  measured  events. 

(e)  the  initial  paramic  estimates  are  not  too  far  from  the 
right  answer.  (Wliat  constitutes  "too  far"  varies  with  the  nature  of 
the  fitting  equations  and  the  measured  data.) 
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*  It  should  be  noted  that  in  applying  Chapman-Kirk  to  the  CD  equations 
of  motion,  Whyte  used  an  unweighted  least  squares  criterion.  Since  the 
angular  and  translational  residuals  of  the  fit  are  not  of  equal  magni¬ 
tude,  Whyte  Das  forced  to  decouple  the  angular  equations  from  the 
translatio^l  equations.  Dissatisfaction  vith  this  enforced  and  often 
unrealistic  decoupting  led  Whyte  and  Hathaway  to  dbundon  an  unweighted 
least  sqi^ea  in  favor  of  a  weighted  maximum  likelihood  criterdon. 

Since  this^  cHterion  was  derived  on  the  assumptiem  of  a  normal  error 
distribution,  their  Maximum  Likelihood  Method  [see  Eeferenaes  35-39  in 
the  Bibliography]  should  yieM  the  'same  final  fit  (albeit  by  a  differ¬ 
ence  path)  as  a  comparably  weighted  least  squares  approach.' 
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,  an  N  by  N  multi- round  matrix. 


the  (j,kl-th  element  of  matrix  A.  [j »k=l, 2, . . .N] 
(a.  .a, ,  a.,  [nondimensional] 

J  J  KK  1 K 

1  +  X  wh'jn  j^k;  a.j^  when  j]i<k 
NR  ^ 

^  ,  an  N-diroeniional  multi-round  vector. 

n«l 


^2  *  system  (4) 

the  k-th  component  of  ?. 

-U 

bj^  [nondimensional] 
the  vector  of  N3  parajneters. 


i  ROMA  input  argument  vector  of  N2  single-round  initial 


c(,t^j^)/e(Pn_i),  FINLIE's  measure  of  convergence, 

the  j-th  parameter.  [j=l, 2, . . .N3] 

<^y  (j»k)-th  influence  coefficient 

evaluated  at  x^,  using  the  current  point  Q. 
[j=l,2,...N2;  k=l,2....N23] 


a  ROME  output  argument  vector,  Eq.  (74). 


a  DUBLIN  output  argument  vector  of  crude  estimates 

*k’  [k=l ,  2, . . . N] 

the  n-th  round  identifier.  [n=l , 2, . . .NR] 

a  DUBLIN  output  argument;  e(P),  the  value  of  e  at  the 
point  currently  stored  in  argument  P. 

a  DUBLIN  input  argument  vector  of  N  adjust-or-hold- 
fixed  flags  associated  with  the  N  paramics  pj^. 

a  COMMON/NAPLES/  input  vector  to  ROME  (and  ROMA)  containing 
the  .N23  single-round  adjust-or-hold-f ixed  flags 
associated  with  the  N23  paramics  qj^. 
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fj  Tj',  system  (1) 

gj  Yy  system  (2) 

h  a  noTidijmensional  positive  constant  in  the  steepest 

descent  technique,  Eq.  (55) . 

IC  the  set  of  multi-round  initial  conditions. 

N  (NR  X  N2)  +  N3,  the  number  of  paramics  in  the  system. 

K'A  N2  X  N23,  the  number  of  influence  equations  for 

systan  (1). 

NB  N1  X  N23,  the  number  of  needed  influence  equations 

for  system  (2) . 

NC  a  DUBLIN  I/O  argument:  originally  zero  (the  "first- 

call"  flag),  it  becomes  the  number  of  paramics  adjusted. 

NF-  a  DUBLIN  input  argument:  1  to  fit  system  (1);  0  to  fit 

system  (2). 

NM  a  DUBLIN  input  argument  vector,  where  NM(J)  is  the 

number  of  data  points  in  the  j-th  round. 

NR  a  DUBLIN  input  argument:  the  number  of  rounds  (hence 

the  number  of  distinct  sets  of  initial  conditions  to  be 
determined) . 

NS  a  DUBLIN  output  argument:  0  means  "CALL  again"; 

1  means  "a  disaster  has  occurred";  2  means  "convergence 
by  FINOiE's  criterion". 

NW  a  DUBLIN  input  argument:  0  to  use  the  user's  weights; 

1  tc  set  all  weights  at  unity. 

N1  a  DUBLIN  input  argument;  the  number  of  measured 

dependent  variables. 

N2  a  DUBLIN  input  argument:  the  number  of  dependent 

variables  in  the  system. 

N3  a  DUBLIN  input  argument:  the  number  of  parameters  in 

the  system. 

N4  the  niunber  of  data  points  fer  all  the  rounds. 
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N23 

P 
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PAR 


paramic 

Pk 


Pk 
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R 

RL 


a  ROME  input  argument:  the  number  of  equations  (N2 
differential  equations  plus  NA  Influence  equations)  to 
be  defined  in  ROME. 

N2  N3,  the  number  of  paramlcs  in  a  single  round. 

(a)  a  set  of  N  multi-round  paramics; 

(b)  a  point  in  the  N-dimensional  paramic  space; 

(c)  a  fJUBTDj  I/O  argiunent  vector  containing  the  N 
paramic  values. 

the  vector  from  the  origin  to  point  P  in  the 
N-dimensional  paramic  space. 

the  value  of  point  P  that  mininiizes  e. 

a  CCMiJON/NAPLES/  input  vector  to  ROME  (and  ROMA) 
containing  the  N3  parameters 

the  value  of  point  P  at  the  end  of  the  n-th  iteration. 
[n=0,l, . . .] 

a  candidate  for  point  P^,  obtained  by  setting 
[n=l,2....] 

parameter  or  initial  condition. 

the  k-th  paramic  (k=l,2,...N)  where  the  order  is: 
first-round  initial  conditions,  second-round  initial 
conditions,...  and  finally  the  N3  parameters. 

(aj^j^)  ^  pj^,  the  k-th  nondimensional  paramic. 

the  single-round  equivalent  of  P. 

the  value  of  point  Q  that  minimizes  y. 

the  k-th  single-round  paramic,  where  the  order  is:  the 
N2  initial  conditions,  then  the  N3  parameters. 

a  DUBLIN  output  argument:  the  N1  by  N4  matrix  of 
residuals,  Eq.  (83),  at  the  current  point  P. 

a  DUBLIN  I/O  argument:  initially  set  at  0.0  to  avoid 
Marquardt's  algorithm;  1.0  to  Invoke  it. 

the  m-th  data  point,  consisting  of  and  the  N1 
dependent  variable  measurements 
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ROMA 

ROME 


RS 


round 


S 

s 

SIG 


s 

u 


u 

w 


X 

XE 

XO 

X 


an  arbitrary  name  for  the  user’s  subroutine  defining 
system  (2) . 

an  arbitrary  name  for  the  user's  subroutine  defining 
system  (ij. 

a  DUBLIN  output  argument  vector  of  N1  error  measures, 
Eq.  (84). 

an  experiment  at  which  measurements  were  taken  and  with 
which  a  distinct  set  of  initial  condition  values  can 
be  associated. 


the  N-dimensional  space  in  which  point  P  has  coordinates 
(Pl»  1?2‘  ’ 


the  N-dimensional  scaled  space  in  which  point  P  has 
coordinates  (pj,  ^2*  Pj^-1 


the  N23-dimensional  single-round  spac6;  in  which  point 
Q  has  coordinates  (qj»q2> • • *1^23^ ' 

a  DUBLIN  output  argument:  the  value  of  s  at  tne 
current  point  P. 

estimated  standard  deviation  of  the  fit,  Eq.  (69). 

the  crude  estimated  standard  deviation  in  paramic  p,  , 
Eqs.  (70-71). 

a  ROME  input  argument  vector,  Eq.  (73);  a  ROMA  output 
argument  vector. 

exp[c^(x-XQ)]  in  system  (4). 

a  DUBLIN  I/O  argument:  the  N1  by  N4  matrix  of  weights 


w,  . 
jm 

the  non-negative  weighting  factor  associatisd  with 


a  DUBLIN  input  argument  vector  of  the  N4  values  x^. 


a  ROME  and  ROMA  input  argument’  the  value  of  x. 


a  DUBLIN  and  ROMA  input  argument:  Xq. 
the  independent  variable  of  the  system 
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the  m-th  value  of  x  at  which  measurements  were  taken. 

the  value  of  the  independent  variable  at  which  all 
initial  conditions  apply. 

(a)  a  vector  of  N2  dependent  variables; 

(b)  a  DUBLIN  input  argument:  the  Nl_by  N4  matrix  of 
dependent  variable  measuraments  y^^. 

a  DUBLIN  output  argument:  the  N2  by  N4  matrix  of 
computed  dependent  variable  values  yj(Xj^,P). 

a  vector  of  N2  single-round  initial  conditions. 

the  j-th  dependent  variable. [j=l, 2, .. .N2] 

the  measured  value  of  y.  at  x 

j  m 

[j=l,2, ...Nl;  ra=l,2, ...N4] 

the  calculated  value  of  at  x^,  using  the  current 
point  P. 

the  N23  by  N23  single-round  matrix  at  point  Q 

the  N  by  N  expansion  of  the  matrix  a  associated  with  round 
En. 


the  (k,n)-th  element  of  a(Q),  Eqs.  (36-38). 

a  vector  point  function  of  Q:  the  vector  in  space  S, 
in  the  direction  of  the  negative  gradient  of  ^  st 
point  Q. 

the  N-dimensional  expansion  of  the  vector  ^  associated 
with  round  En. 


Ip 

n 

y(Q) 


the  k-th  component  of  [k=l , 2, . . .N23] 

the  increment  vector  in  S  from  point  P  to  the  nearby 
point  P^^j. 

the  increment  vector  in  S,  from  point  Q  to  the  nearby 
point 

the  single-round  equivalent  of  e(P),  Eq.  (14). 


LIST  OF  SYMBOLS  (continued) 

e{F)  the  nondimensional  sum  of  the  weighted  squares  of  the 

residuals,  Eq.  (7);  the  value  of  the  scalar  point 
function  e  at  point  P. 

X  a  nonnegative  constant  added  to  each  diagonal  element 

of  the  scaled  form  of  matrix  A  and  adjusted  by 
Marquardt's  algorithm  so  that  <  e(P^). 

XnA»  consecutive  trlAl  values  assigned  to  X  in  an  effort 

to  move  from  point  P^^,  where 


Superscripts 

C  )  a  row  vector 


(  ) 
( 


the  transpose  of  a  row  vector;  that  is,  a  column  vector 
the  scaled  (hence  nondimensional)  form  of  (  ) . 
d(  )/dx 


Subscripts 

t  Id 

(  )c 


(  ); 
(  ) 


Si 


denotes  the  dimensions  of  [  ] 
the  components  of  (  )  are  in  space  S 

the  components  of  (  )  are  in  space  S 

the  components  of  (  )  are  in  space  Sj. 
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f-(J>  ■  l.o  IF  THE  CURRENT  VALUE  OF  P(J)  IS  TO  BE 
adjusted  iV  THE  FITTING  PROCESS 
NV  a  THE  WEIDHT  FLAG  ASSOCIATED  WITH  W  BELOWt  WHERE 

NW  a  0  IF  the  USER'S  WEIGHTS  IN  ARGUMENT  W  ARE  TO 
BE  USED 

NW  a  I  IF  all  WEIGHTS  ARE  UNITY,  IN  THIS  EVENT*  THE 
THE  FIRST  TIME  THIS  SUBROUTINE  IS  CALLED.  ALL 
elements  of  matrix  W  are  set  TO  1,0.  (HENCE. 

the  user  need  not  establish  w  when  all  its 
elements  are  l.o.) 

w  a  THE  Nl  BY  N4  MATRIX  OF  WEIGHTS  ASSOCIATED  WITH  INPUT 
Y  ABOVE,  SEE  ARGUMENT  NW  ABOVE. 

INPUT/OUTPUT  ARGUMENTS  ... 
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. . PINE.'  «  I.C.  FOR  FIRST  ROUND 

. . .  a  I.C.  FOR  SECOND  ROUND 

P(NR*N2-N2''1)  ....P(NR‘*N2)  a  I.C,  FOR  LAST  ROUND 
P(NR*Nc-l> ,.,,,,,P{N)  a  PARAMETERS 
RL  a  A  MARQUARCT  ARGUMENT.  BEFORE  DUBLIN  IS  CALLED  THE 
FIRST  TIME. 

SET  RL  a  0*0  IF  THE  MARQUARDT  ALGORITHM  IS  TO  BE 
OFITTEO.  THEREAFTER.  RL  WILL  REMAIN 
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SET  RL  a  1.0  IF  THE  MARQUARDT  ALGORITHM  IS  TO  BE 
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THIS  SUBROUTINE  THEN  ESTABLISHES  NC  AS  THE  ACTUAL 
NUMBER  OF  INITIAL  CONDITIONS  AND  PARAMETERS  BEING 
determined  (1  ,L£.  NC  .LE.  N> 
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R(I.J)  a  Yi  I.j)  -  YC(I.U) 
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SUM  OF  the  SQUARES  OF  THE  RESIDUALS  OVER  ALL  THE 
POINTS.  OVER  ALL  THE  MEASURED  DEPENDENT  VARIABLES 
AND  OVER  all  THE  ROUNDS. 

SIG  a  THE  estimated  STANDARD  DEVIATION  OF  THE  FIT 
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NS  a  0  IF  The  process  has  not  YET  CONVERGED  BY 

ThE  CRITERION  BUILT  INTO  SUBROUTINE  DUBLIN 
a  I  IF  all  OUTPUT  ARGUMENTS  ARE  INVALID  (PROBABLY 
BECAUSE  SOME  INPUT  ARGUMENTS  ARE  TOO  LARGE 
FOR  CERTAIN  DIMENSIONED  ARRAYS),  THE  CALLING 
PROGRAM  SHOULD  TAXE  SPECIAL  ACTION  (E,3.. 
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FORM  DPM  ■  THE  SQUARE  OF  THE  MAQiMITUDE  OF  THE 
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EXPAND  PC  TO  FULL  SIZE  AS  PB.  THE  CANDIDATE 
replacement  for  input  point  P. 
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PB<J)  *  P{J) 


160 

170  CONTINUE 
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EB  (AND  associated  ARRAYS  YCt  Rt  RS*  ALPHA  AND  BETA) 
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CALL  LONDON  ( POME  * NF » M 1 , M2 tMa* M7 »M8» 60 f MR, NM , X 0 t X . Y » W ♦ E , PB . 
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IF  THE  ANGLE  BETWEEN  BATA  AND  DELTA  P  IS  LESS  THAN 
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length  OF  delta  P  and  GO  HACK  TO  PART  4, 
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no  190  J  ■  i.LO 
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POO  CONTINUE 
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P(J)  ♦  O.I»(PB(J)-P(J)) 


increase  MARQUARDT'S  lambda  and  GO  BACK  TO  PART  3, 
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THE  HAROUARDT  ITERATIVE  PROCESS  HAS  SEEN  COMPLETED 
satisfactorily,  update  error  measures  EPS  AND  SIG 
TEST  FOR  CONVERSENCE,  UPDATE  POINT  P  AND  COMPUTE 
ERROR  estimates  EK. 


C 

c 

c 


c 

c 

c 

L 

c 

c 

c 
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210  CONTINUE 
RL  »  QL 
EPS  »  EB 

SIG  »  SQRT(EB/EM) 

CR  B  1.0  -  EB/EA 

IF  (CR  .GE.  0.  .AND.  CR  .LT.  0,000010)  NS  ■  2 

K  B  1 

DO  220  J  ■  1,LD 
P(J)  »  PB(J) 

IF  (F(J)  .EQ.  0.)  GOTO  215 

EK(J)  B  SIG*S(K)*SQRT(QAMMA(K,K)*D1AG) 

K  B  K  ♦  1 
GOTO  220 

21S  EK(J)  B  0. 

220  CONTINUE 
RETURN 
END 


subroutine  LONDON  ( ROME , NA ♦ N I » N2, N3, N? t N8 * NL ♦ NR, NM , X A, X , Y ♦ V, F , P» 
1  EPS, NC. ALPHA, BETA, YC*R*RS, IS) 


»  « 


«  «  *  * 


* 

« 


FOP  A  GIVEN  SET  OF  PARAMETER  AND  I,C.  VALUES  AND  A  GIVEN 
MULTI-ROUND  SET  OF  MEASUREMENTS#  THIS  SUBROUTINE  PRODUCES 
THE  ERROR  MEASURE  E^S,  THE  COMPUTED  DEPENDENT  VARIABLE 
VALUES,  the  RtSlUUALS  AND  THE  ALPHA  AND  riETA  ARRAYS  POR  THE 
multi-round  OaTa,  AL>  ARGUMENTS  ARE  DEFINED  IN  THE  COMMENTS 
VITHIN  SUBROUTINE  DUBLIN, 


DIMENSION 


1 


dimension 
the  above  or 

DIMENSION 

the  above  01 

DIMENSION 
THE  ABOVE  01 
DIMENSION 
THE  ABOVE  01 
DIMENSION 

tmk  above  DI 

DIMENSION 
the  ABOVE  01 

dimension 
the  above  01 


NM(  1)  ,P(  1)  ,F  {  n  ,X<  1)  ,Y(N7,n  ,V(N7,  1)  . 
YC(N««1)  ,R(N7,1)  ,RS(1)  .ALPHA  (NLtl!  tSETAd) 
RSQ(IO) 

MENSION  ASSUMES  THAT 
C(20) ,CE (20) 

MENSIONS  ASSUME  THAT 
CP(40) ,FP(40) 

MENSIONS  ASSUME  THAT 
ALFA(3600> .BATA (60) 

MENSIONS  ASSUME  THAT 
RE(200) ,VR(200) ,YM(200) 

MENSIONS  ASSUME  THAT  Ni 
YCOMP ( *00 ) 

MENSION  ASSUMES 
XX (1000) 

MENSION  ASSUMES 


Nl  .LE,  10 
N2  ,LE.  20 
N3  ,LE.  *0 
N2  ♦  N3  .LE.  60 

•(MAX,  ELEMENT  OF  NH)  ,LE,  200 

That  n2*(max,  element  of  nm)  ,le,  *oo 


external  romp 
e**  pftRT  1.  *•»  preliminaries 


THAT 
NM  (  1  ) 


♦  NM(2)  ♦ 


♦  NM(NR)  .LE.  1000 


FNLE 

?99 

FNLE 

300 

• 

fnle 

301 

fnle 
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fnle 

303 

FNLE 

30* 

FNLE 

305 

FNLE 
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FNLE 

307 

FNLE 

30H 

FNLE 

309 

FNLE 
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FNLE 

311 

FNLE 

312 

ENLE 

313 

fnle 

31* 

fnle 

315 

ENLE 

316 

fnle 

317 

ENLE 

31H 

ENLE 

319 

ENLE 

320 

ENLE 

321 

- 

-  - 

32? 

ENLE 

3?3 

'I 

ENLE 

32* 

ENLE 

325 

fnle 

326 

* 

enle 

327 

ENLE 

328 

# 

enle 

329 

» 

enle 

330 

enl: 

331 

• 

enle 

3  32 

# 

ENLE 

333 

# 

enle 

33* 

• 

enle 

335 

• 

enle 

3  36 

ENLE 

337 

enle 

33B 

enle 

3  39 

ENLE 

3*0 

enle 

3*1 

enle 

3*? 

enle 

3*3 

ENuE 

3** 

enle 

3*5 

enle 

3*6 

enle 

3*7 

ENLE 

3*fl 

enle 

3*9 

enle 

350 

enle 

351 

enle 

35? 

enle 

35  3 

enle 

enle 

35* 

3S5 

enle 

3  56 

enle 

357 

enle 

3  5  A 

82 


o  o  o 


'^10 
?2  0 

230 

24  C 

««» 


250 

250 


270 

1 


IS  s  0 

FNLE  355 

MR  B  NR 

F.NLE  360 

Ml  s  N 1 

FNLE  361 

M2  «  N2 

FNLE  362 

M3  ■  N3 

FNLE  363 

M23  »  M2  «  M3 

FNLE  364 

M6  ■  M2*MR 

FNLE  365 

LC  a  1  +  M6 

FNLE  366 

LD  «  M3  ♦  MS 

FNLE  367 

JX  «  1-(14.M23)*(LC-M2) 

FNLE  36ft 

no  220  N  ■  1»LD 

FNLE  369 

ALPHA(NiN)  >  0. 

FNLE  370 

BETA(N)  ■  0, 

FNLE  371 

IF  <N  .EQ.  LD)  goto  220 

FNlE  37? 

MA  s  N  ♦  1 

FNLE  373 

DO  210  K  «  MA«LO 

FNuE  374 

ALPHA(K.N)  r  0. 

FNLE  375 

ALPHA(N.K)  ■  0, 

FNLE  376 

CONTINUE 

FNLE  377 

CONTINUE 

FNLE  378 

DO  230  I  ■  1,M1 

FNLE  379 

RS(I)  ■  0. 

FNLE  580 

CONTINUE 

FNLE  381 

MC  =  0 

FNLE  382 

DO  240  K  =  1,M3 

FNLE  383 

KA  =  K  ♦  M6 

FNLE  384 

CP(K)  a  R(KA) 

FNLE  385 

FP(K)  s  F(KA) 

FNLE  386 

IF  (FP(K)  ,NE.  0.)  MC  a  MC  ♦  1 

FNLE  387 

CONTINUE 

FNLE  388 

JA  a  0 

FNLE  l-i89 

JH  S  0 

FNLE  390 

EP  *  0. 

FNLE  391 
FNLE  39? 

PART  2.  ***  THE  DO-LOOP  FOR  HANOLINO  MULTIPLE  ROUNDS 

FNLE  393 
FNLE  394 

no  370  JR  a  i,MR 

FNLE  395 

M4  a  NM(JR) 

fnLF  396 

DO  260  M  a  l,M4 

FNLE  397 

lA  a  (M-1)<*M1 

FNLE  398 

JA  a  JA  ♦  1 

FNLE  394 

XX  (M)  a  X  1  JA) 

4NLE  40n 

DO  250  I  a  1,M1 

FNLE  401 

IM  a  lA  +  I 

FNLE  40? 

YM ( IM)  a  Y ( I ,JA) 

FNLE  403 

MR  (IM)  a  Md.JA) 

NL  E  4  0  4 

CONTINUE 

'■NlE  4  05 

CON'^INUE 

FNLE  406 

LA  a  ( JR-l ) *M2 

FNLE  4(^7 

no  270  K  a  1 ,M2 

FNLE  408 

LA  a  la  ♦  1 

FNLE  409 

C(K)  a  P(LA) 

FNLE  410 

CF  (K)  a  F  (LA) 

FNLE  411 

IF  (CF  (K)  ,NE.  0.)  MC  a  MC  ♦  1 

F  N  L  E  4  1? 

CONTINUE 

FNlE  413 

CALL  PARIS  (ROME.NA, lfMl,M2,M23,M4,C.CF,CPtFP,MR,XA.XX,YMv 

FNLE  4  1  4 

YCOMP.REtRSQ. ALFAtBATAt IR) 

FNLE  415 

IF  (IR  ,F(3.  0)  GOTO  2fl0 

FNLE  416 

IS  a  1 

FNLE  4  1  ; 

PRINT  275 

FNLF-  418 

83 


W!*IPWJI^"||." 


'•/5 


2ao 


FT>«><ATaH  aOi^.'UNSUCCESSrUL  WETuaw  FROM 
I'/IH  .lOX.'ALL  subroutine  DUBLIN  OUTPUTS  INVALID. 

return 

CONTINUE 
DO  200  1  »  l*Nl 

RS(n  «  HSSli  ♦  RSQd) 


SUBROUTINE  PARIS 
* ) 


2'iO 


300 

310 


Z20 


330 


34  0 


35  0 


EP 

CONTINUE 
DO  3X0  M 
NI  • 

NJ  » 

J3  » 

00  300  J 
IM  B 


EP  •i'  RS8(I) 


»  I  *NA 
(W  1  It  *Ml 
(M-l»  *M2 
JB  ♦  1 

■  1»M2 
NI  ♦  J 


JM  ■  NJ  ♦ 

YC(J*JB)  « 
IF  (J  .LE. 
CONTINUE 

continue 

LA  ■  1  ♦  (.)R-n«Mt 
LB  ■  LA  ♦  M2  1 
J  «  1 
JJ  «  1 


YCOMP'.UMl 
HI)  R(JfJB) 


RE(IM) 


DO  340  N 


.A. LB 


If  «  N 

J  «  J  ♦  K 

-  LA 

ALPHA (N»KI 

»  alfa(J) 

IF  (K  ,6T. 

N)  A..PHA(K.N) 

J  ■  J  ♦  1 

K  ■  K  ♦  1 
IF  (K  .LE. 

L8)  goto  320 

K  a  LC 
ALPHA (N»K) 

«  ALFA(J) 

ALPHA (K, N) 

■  ALFA(J) 

J  ■  J  ♦  1 

K  a  K  •  1 
IF  (K  ,LE. 

LU)  goto  330 

BETA(N)  » 

BaTA ( JJ) 

JJ  ■  JJ  ♦ 

1 

ALFA( J) 


continue 

JJ  *  1  ♦  M? 

DO  360  N  B  ‘.CfLD 
K  »  N 

J  a  JX  ♦  (1 

ALPHA (N*K) 

IF  (K  ,ST.  N  tANO,  JR  .EQ»  HR)  ALPHA(K»N)  ■  ALPHA(N»K) 
J  *  J  ♦  1 
K  »  If  ♦  1 
IF  (K  ,LE.  LD) 

BETA(N)  =  RETAIN) 


♦  Mp3)*N 

*  Alpha (NfK)  ♦  alfa(j) 


GOTO  350 

BATA ( JJ) 


J  J  «  J  J 

360  CONTINUE 

370  continue 
EPS  ■  EP 
NC  a  MC 

return 

END 
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4  1  N 
4?0 
421 
42? 
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4  26 
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4  36 
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4  3fl 
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SUHOOUTINE  PARIS(ROM£..NA,NBtNl.N2*N23»N<^(,C.CE.‘>»>‘-'E.«»X<k.X,  YM. 
1  YC.R»RSQiALEA,BATAt IRJ 


*  FOR  A  GIVEN  S;iT  OF  PARAMETER  AND  IC  ESTIMATES  AND  TH£  GIVEN  * 

*  MEASUREMENTS  FOR  A  SINGLE  ROUND#  THIS  SUBROUTINE  PRODUCES  * 

*  THE  COMPUTED  DEPENDENT  VARIABLE  VALUES#  THE  RESIDUALS  AND  • 

*  IF  NB  .NE,  0)  THE  ALPHA  AND  BETA  ARRAYS  FOR  THE  GIVEN  ROUND.* 
«  * 


# 

* 

# 

» 


# 


# 

« 

« 


# 

# 

# 


* 


INPUT  ARGUMENTS  ... 

POME  •  THE  dummy  name  OF  ThE  SUBROUTINE  {WRITTEN  BY  THE 
USER)  that  DEFINES  EITHER 

{A)  THE  dependent  VARIABLES  (NA  »  0) 

OR  <B)  THE  DERIVATIVES  OF  THE  DEPENDENT  VARIABLES 
(NA  .NE.  0) 

NA  ■  THE  FLAG  THAT  INDICATES  THE  NATURE  OF  SUBROUTINE 
ROME  (SEE  PREVIOUS  ARGUMENT) 

NB  ■  0  IF  arguments  ALFA  AND  BaTA  BElOB  ARE  NOT  TO  BE 
COMPUTED 

■  1  IF  ALFA  AND  BATA  ARE  TO  BE  COMPUTED 
•  mUMBER  of  measured  dependent  variables  (Nl.GE.l) 

N?  «  TOTAL  NUMBER  OF  DEPENDENT  VARIABLES  (N2.GE,N1> 

N23  •  total  number  of  parameters  and  initial  conditions 

(N23.GE.N2) 

N4  s  number  of  measurements  taken  on  each  of  the  N1 
MEASURED  variables 

C  ■  VECTOR  OF  The  N2  INITIAL  CONDITION  ESTIMATES 

CF  «  VECTOR  OF  T ME  N2  INITIAL  CONDITION  ^^LAGS 

(0.  IF  the  INITIAL  condition  VAlJE  IS  FIXED) 

P  ■  VECTOR  OF  The  N3  parameter  VALUES 

PF  c  VECTOR  OF  ThE  N3  PARAMETER  FLAGS 

(0.  IF  the  parameter  VALUE  IS  FIXED) 

M  ■  VECTOR  OF  TmE  M  BY  N4  MEIOmTS  ASSOCIATED  HITH  INPUT 
YM  defined  8EL0M 

XA  «  THE  reference  X  VALUE  AT  <HICH  INITIAL  CONDITIONS 
apply 

X  s  vector  of  The  N4  values  of  the  INOERENDENT  variable 
AT  RHICH  HEaSUREMENTS  MERE  TAKEN 
YH  >  VECTOR  OF  TmE  N1  BY  N4  MEASURED  Y  VALUES  FOR  ONE 
ROUND.  MhERE 

YM(IM)  ■  measured  VAlUE  of  Y(I>  AT  X(M) 

IM  ■  1  .  (R-I)*N1 

I  -  1.2.  ...  N1 
M  >  1.2,  ...  N* 


OUTPUT  arguments  ... 

YC  ■  VFCTOS  OF  N2  BV  N4  COMPUTED  v  VALUES  FOR  ONE  ROUND, 


R  s 


RSQ 


MHFRE 

YC(JM)  «  COMPUTED  VALUE  OF  Ylj)  aT  X(M) 

JM  ■  J  ♦ 

J  •  1.2.  ...  N2 

VECTOR  OF  Nj  BY  NA  RcSIOUAlS.  MhERE 
R(  IM)  a  VMdM)  -  YC(  IM2)  . 

IMJ  •  I  ♦  (M-1)»N2 

vector  of  N1  error  measures,  where 

»son)  ■  wlighteo  sum  of  thf  squares  of  the  residuals 

!N  The  I-Th  DEPENDENT  VARIABLE 
VECTOR  OF  N23  BY  N23  ALPHA  VALUES#  MH-RE 
ALF4(Nr)  »  AlPMA(N,K) 
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*  NK  “  N  ♦  <K-U*N23 

*  BATA  a  VECTOR  OF  N23  BETA  VALUES 

*  IR  »  I  IF  A  warning  was  printed  ANO  A  RETURN  EXECUTED 

*  BEFORE  03TAINING  THE  ABOVE  OUTPUT 

*  *0  OTHERWISE 

* 


DIMENSION  Cn»  tCF(  1  !  .P(l )  ,PF{1)  .K(l»  .yM(  1)  ♦TC(  1»  ♦«  (I )  .«S0(  U  * 
1  ALFA(U  tBATAd)  iWd) 

DIMENSION  U(400)  (UU(AOO)  tSdCOOO) 


•  *«••*•#••••***•••*****•••*••••*• 

« 

*  WARNING  ...  THE  DIMENSIONS  OF  U  ANO  DU  ABOVE  MUST  EQUAL 

*  OR  exceed  N5  ■  N2*d^N23»»  THE  NUMBER  OF  EQUATIONS 

•  IN  subroutine  ROHE.  THE  DIMENSION  OF  ARRAY  S 

*  ABOVE  MUST  EQUAL  OR  EXCEED  N4»N5.  FOR  EACH  VALUE 

•  OF  the  independent  VARIABLE  X, 

• 

*  U  subscripts...  DENOTE  THE  COMPUTED  VALUES  OF  ... 

• 

•  1*2....  N2  THE  N2  DEPENDENT  VARIABLES 

•  N2*l«.,.  N5  the  N2*N23  PARTIAL  DERIVATIVES.  WHERE 

•  PARTIAL(U«<»  »  UIU*K*N2) 

« 

•  ARRAY  DU  denotes  DU/OX.  ARRAY  S  IS  A  COMPILATION 

•  OF  all  The  u  values,  that  is, 

•  Sd).,..S(N5)  •  Ud)  ,„..U(NS)  at  X  d) 

•  S(N5*l)  ,..,S(2PN5)  •  Ud),...J(NS)  AT  X(2),  ETC. 

•i 


common  /NAPLES/  PAR(<.0)  .FLAG(60) 


•  TmF  ABOVE  labelled  COMMON  MyST  BE  LINKED  TO  ANO  USED  IN  THE 

•  USER'S  subroutine  ROME.  THE  VALUES  ARE  ESTABLISHED  HERE'  IN 

•  SUBROUTINE  PARIS. 

•  par  ■  The  ARRAY  OF  N3  PARAMETER  VALUES 

»  Flag  «  the  array  of  n23  flags  for  a  given  rouho,  where 

•  FlA6(K)  *  K-TH  initial,  condition  flag,  K  ■  1.2...N2 

•  FLAG(K.N2)  »  K-Tm  parameter  flag,  K  ■  1,2, ..N3 


external  ROME 

data  HMIN.O.OELX/- 1 , 0,  .00010,  ,12S/ 

•••  part  d)  preliminaries 
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N4 
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M?* ( 1 »M?3 ) 

IF 
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IR  »  1 
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S66 
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ft 

fnle 

56H 
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fNlE 

569 

ft 

fnle 

570 

FNlE 

S71 

"NLE 

572 

fnlE 

5  73 

ft 

fnle 

574 

ft 

fnle 

5Y'-, 

ft 

fnle 

576 

ft 

fnle 

577 

ft 

fnle 

57H 

ft 

fnle 

‘•79 

c 

fnle 

5H0 

ft 

fnle 

5P1 

ft 

fnlE 

5h2 

ft 

fnle 

5fi3 

ft 

fnle 

554 

fnle 

505 

fnle 

506 

fnle 

507 

fnle 

500 

fnle 

509 

F'JLE 
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fnle 
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fnle 

592 

fnle 
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fnle 

594 

fnle 

59S 

fnle 

S9b 

fnle 

597 

fnle 

59m 
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n  n  ry  n  nnoonoo 


PRINT  10 

fnlE 

s  'yn 

10 

FORMATdHO.*  SUBROUTINE  PARIS  WARNING  ») 

FNLE 

■soo 

PRINT  22. M5 

FNLE 

601 

FORMAT  IIH  »10X,  •  INCREASE  DIMENSIONS  OF  !J  AND  DU  (ANO  IN  SU8R0UFNLE 

602 

ITINF  MERSO.  INCREASE  DIMENSIONS  OF  T.  G  AND  S)  TO  *.I5) 

FNLE 

6  03 

2* 

Mfe  .  M4*M5 

FNLE 

6  0  4 

IF  (M6  ,LE.  10000)  GOTO  28 

FNLE 

60S 

IR  ■  1 

FNLE 

60  6 

PRINT  10 

FNLE 

607 

PRINT  26. Mi 

FNLE 

60B 

2b 

FORMATdH  .lOX.  »  INCREASE  DIMENSION  OF  S  TO  *.I5) 

FNLE 

609 

2a 

CONTINUE 

FNLE 

6  i  0 

IF  (IR  .EQ.  1)  RETURN 

FNLE 

611 

IH  «  0 

FNLt 

612 

H  «  DELX*(X(M4)-Ail)  )/FL.OAT(M4-l) 

FNLE 

6  11 

JA  B  M?  ♦  1 

FNLE 

6  1  4 

no  30  JB  a  JA.MS 

fnlf 

61-) 

U(JB)  ■  0,0 

FNLE 

6  1  6 

DU(J8)  >  0.0 

FNLE 

617 

30 

continue 

FNLE 

61fl 

no  32  K  m  1.M2 

'=‘NLE 

6  1  ■» 

FLAG(K)  >  CF(K) 

FNLE 

620 

32 

CONTINUE 

FNLE 

621 

M3  ■  M?3  -  M2 

enle 

6?? 

no  34  K  «  I.M3 

enle 

623 

PAR(X)  a  P(K) 

enle 

624 

KA  ■  K  .  M2 

FNLC 

62S 

FLA6{KA)  a  PF(K) 

FNLE 

626 

34 

CONTINUE 

enlE 

FNLE 

627 

626 

PART  (2)  DETERMINE  KU  •  NO,  OF  POINTS  IN  X  LESS  THAN  XA 

FNLE 

629 

AND  SR  a  NO,  OF  POINTS  IN  X  GREATER  THAN  XA, 

FNLE 

6  10 

IF  XA  COINCIDES  WITH  A  POINT  IN  X.  THEN  THE  'COMPUTED* 

FNLE 

631 

Y  VALUES  AT  THAT  POINT  ARE  THE  INITIAL  CONDITIONS  FROM 

FNLE 

63-^ 

SUBROUTINE  BONN, 

FNLt 

fnlE 

6  33 

6  3  4 

no  40  M  s  1 , M4 

FNLE 

4  3*^ 

IF  (XA-X(M))  70.50,40 

FNLE 

'  14 

40 

CONTINUE 

fnlE 

'  37 

KL  a  M4 

FNLE 

63H 

KR  a  0 

FNLE 

6  39 

GOTO  80 

fnlE 

6*.  ) 

SO 

KL  a  M  -  1 

FNLE 

64  1 

KR  a  M4  -  M 

fnlE 

64? 

IM  a  1 

FNLE 

6  4  3 

call  B0NN(H2,MS,C. Jl 

fnlE 

6  4  <* 

LB  a  KL^MS 

FNLE 

64S 

no  60  L  a  1,MS 

FNLE 

6  46 

LB  a  lb  ♦  1 

FNLE 

647 

S(L3)  a  U(L) 

FNLE 

6  4  <3 

60 

CONTINUE 

FNLE 

649 

GOTO  80 

enle 

6S0 

7C 

KL  a  K  -  1 

enlZ 

6  S  1 

KR  «  M4  -  KL 

fnle 

enle 

6  62 

6S3 

•  •  • 

PART  ii',  FOR  EACH  POINT  IN  X,  SOlVE  THE  SYSTEM  OF  EQUATIONS 

fm.  £ 

6S4 

DEFINED  IN  SUBROUTINE  ROME  TO  OBTAI-!  COMPUTED  Y  VALUES. 

'  N  L  E 

6SS 

TmESE  Y  values  are  STORED  IN  S. 

FN^E 

6  6  S 

ao 

CONTI NUE 

FNwE 

bS  ’ 

IE  (KL  .fQ.  0)  GOTO  90 

FNL*- 

66« 

o  o  o  non 


Hl’r«-irtv<rtn<5^»«  >fr^"»i*i/«Fififtl'*'r'"“<"ii»/:xv''  *»)  r"' 


''™'"""|.  iimwimiiiimiwiipi 


no 


I A  »  -1 

FNLE 

6S9 

IH  *  0 

FNLE 

660 

LK  =  KL 

FNLE 

661 

ROTO  ICO 

FNLE 

66? 

lA  »  1 

FNLE 

663 

IH  a  M4  1 

FNLE 

664 

LK  »  IH  -  KR 

FNLE 

66S 

XI  »  XA 

FNLE 

666 

IF  (IM  ,EQ.  0)  CALL  BONN (M2.M5.C.U) 

FNLE 

667 

IM  o  0 

FNLE 

668 

a  X!LK) 

FNLE 

6bR 

ir  <NA  .NE,  0)  SOTO  114 

FNLE 

670 

CALL  R0ME(C.X1'.X2.U) 

FNLE 

671 

Q(UO  116 

FNLE 

67? 

CONTINUE 

FNLE 

673 

HA  =  0ELX*ARS(X2-X1) 

FNLE 

674 

HR  X  ABS(HJ 

FNLE 

67S 

H  *  AMINKHA.HR) 

FNLE 

676 

CALL  MERSO  (R0ME.M5»Rl.X2«Ue0U.H.HMIN.(J) 

FNLE 

677 

CONTINUE 

FNLE 

678 

LB  =  (LK-1)»M5 

FNLE 

67R 

DO  l?n  L  X  l.MS 

FNLE 

680 

LB  *  LB  ♦  1 

fnle 

6B1 

S(LB)  *  U(L) 

FNLE 

68? 

CONTINUE 

fnle 

683 

LK  a  LK  ♦  I A 

fnle 

664 

IF  (LK  .NE,  IB)  GOTO  HO 

fnle 

68S 

IF  (lA  .EQ.  1)  GOTO  130 

fnle 

6  86 

IF  (KR  .NE,  0)  GOTO  RO 

fnle 

FNLE 

687 

688 

*»»  PART  (4) 


CONVERT  VECTOR  S  TO  VECTOR  YC  AND  FORM  THE  RESIDUAL 
VECTOR  R. 


130  CONTINUE 
DO  ISO  M 
NR  « 


«  nM4 
(H-l ) •Ml 
NY  i  (M-1)*M2 
LA  ®  (M-1)*M5 
00  140  J  *  i»M2 

IM  »  NR  ♦  J 

JM  «  NY  ♦  J 

LR  s  LA  ♦  J 

YC(J«)  =  S(L8) 
IF  (J  .LE.  Ml 
IF  (J  ,Lt.  Ml 
140  CONTINUE 

ISO  CONTINUE 


.AND.  M(IM) 
.AND.  M(IM) 


.NE. 

.EO. 


0.  ) 

0.) 


RilM) 
R  (I»«) 


YM( TM) 

0. 


•  **  part  (S’/ 


DO  170  I 


COMPUTE  VECTOR  RSQ 


1  .Ml 
RM  s  0.0 

no  \40  M  a  1.M4 

IM  -  I  ♦  (M-l)*Ml 

RM  «  RM  ♦  W ( IM) •R ( IM> ••£ 


FNLE  68R 
FNLE  bRO 
FNLE  (SRI 
FNLE  6R2 
FNLE  6R3 
FNLE  bR4 
FNLE  6RS 
FNLE  SRb 
FNLE  6R7 
FNLE  bRB 
FNLE  SRR 
FNLE  700 
YC(JM)ENLE  701 
FNLE  70? 
FNLE  703 
FNLE  704 
FNLE  70S 
FNLE  70b 
FNLE  7C7 
FNLE  708 
FNLE  TOR 
FNLE  710 
FNLE  711 
FNLE  71? 


160 

CONT  INiJE 

fnle 

713 

RSQ( n  »  RM 

fnle 

714 

170 

CONTINUE 

fnle 

7  IS 

IF  iN8  ,EO.  0)  REfURN 

FNLE 

7  1  b 

r 

FN!  ' 

717 

C 

•  »  :> 

PART  (6)  compute  VECTOR  ALFA 

fnlE 

7  18 

88 


^*<»<^oonooooooo  ooo 


fnle 

719 

DO  210  K  «  l*M23 

fnle 

7?0 

LA  ■  K*M2 

FNLE 

721 

NL  ■  <K-1)*M23 

FNLE 

72? 

00  200  N  ■  K«M23 

fnle 

7  23 

LB  ■  N»M2 

FNLE 

724 

NK  ■  NL  ♦  N 

fnle 

72S 

ALrA(NK)  ■  0.0 

fnle 

7?h 

DO  190  I  ■  l*Ml 

FNLE 

727 

ALF  »  0.0 

fnle 

72P 

LC  ■  LA  ♦  1 

fnle 

729 

LD  »  LB  ♦  I 

FNLE 

730 

00  180  M  >  1 fHA 

fnle 

731 

IH  ■  I  ♦  (M-U*M1 

fnle 

7  3? 

ALF  ■  ALF  ♦  W( m)*S(LC)*S(LD) 

fnle 

7  33 

LC  ■  LC  •*  MS 

fnle 

T34 

LO  ■  LO  ♦  MS 

fnle 

73S 

180 

CONTINUE 

FNLE 

73b 

ALFA(NK)  at  ALFA(NK)  ♦  ALF 

fnle 

737 

190 

CONTINUE 

fnle 

73H 

IF  (N  .EQ.  K)  GOTO  200 

fnle 

739 

KN  ■  K  ♦  (N-1)*M23 

fnle 

740 

ALFA(KN)  ■  ALFAINK) 

fnle 

741 

200 

CONTINUE 

fnle 

74? 

210 

CONTINUE 

fnle 

74  3 

fnle 

744 

•  «» 

PART  (7)  COMPUTE  VECTOR  BATA 

fnle 

74S 

fnle 

74b 

DO  240  N  r  1,M23 

fnle 

747 

LA  »  N*H2 

fnle 

74« 

8ATA(N)  >  0.0 

fnle 

749 

00  230  I  ■  l.Ml 

fnle 

7S0 

BAT  ■  0.0 

fnle 

7S1 

!.B  ■  LA  ♦  I 

fnle 

7S? 

DP  220  M  »  lfM4 

FNLE 

7^3 

IN  »  I  ♦  (M-l)»Ml 

fnle 

7S4 

BAT  «  BaT  ♦  W(IM)*R(IM)*S{LB) 

fnle 

7SS 

LB  ™  L8  ♦  MS 

fnle 

7Sb 

220 

CONTINUE 

fnle 

7S7 

BATA(N)  >  RATA(N)  ♦  BAT 

FNLE 

7bB 

230 

CONTINUE 

fnle 

7b9 

240 

CONTINUE 

fnle 

7  (SO 

RETURN 

fnle 

7bl 

END 

fnle 

7b? 

fnle 

7b  ^ 

-  _ 

- 

-  - 

7b4 

fnle 

7bb 

SUBROUTINE  BONN ( N2 , NS . C . U ) 

fnle 

7hb 

fNlE 

767 

• 

fnle 

7bB 

» 

• 

fnle 

7  69 

• 

THIS  SUBROUTINE  (CALLED  BY  SUBROUTINE  PARIS)  ASSIONS  INITIAL 

• 

f  t  L  f 

■'TO 

• 

CONDITIONS  (THE  VALUES  AT  X  •  XA)  TO  VECTOR  U. 

FOR  THE 

% 

fnl' 

7  7  ■ 

• 

DEFINITIONS  OF  THE  ARGUMENTS.  SEE  THE  COMMENTS 

IN  subroutine 

FNLt 

y 

• 

PARIS. 

♦ 

FNLE 

7  1 

# 

fnlE 

774 

• 

• 

fn.  f 

77b 

FNLt 

776 

dimension  C(N2) .U(NS) 

PNLt 

7  7  7 

M2  •  N2 

fnle 

7  7h 
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rt  nn 


MS  B  N5 

no  10  j 

U(  J1 

10  CONTINUE 
LA  «  0 


DO  TO  JA 
LA  « 


•  l»N2 
LA  ♦  M2 

DO  20  JB  «  1»M2 
LB  »  LA  ♦  JB 
IE  (JA  ,EQ.  JB)  U(L8) 
IF  (JA  ,NE,  J8)  U{L8) 
CONTINUE 
30  CONTINUE 

IF  (LB  .QEt  MS)  goto  50 
LA  =  LB  ♦  1 
no  40  LB  a  LA»M5 
U(LB)  «  0,0 

40  continue 
SO  continue 

RETURN 

END 


1.0 

0.0 


SUPROUTINE  MFRSO  ( FUNC , N , X , Z . Y , F , H , HMI N. E ) 
dimension  Y(l),F(l)  $  DIMENSION  T ( 40 0) t 6 (♦ 00 ) , S ( 40 0 ) 

LOGICAL  BC.8E.flH.8R.eX  S  NTaNS  ZTaZS  HMIaHMiNt  ET»ABS(E) 
IF(HMl.LT.0.0)HMlaO,0l«ABS(H)  %  BH*BRaBX« . TRUE . 
flCaO.O.LT.ET.AND.ET. LT.l.O  $  E5aET*5.0 

IF  (  (ZT.GT.X.ANO.H.LT.O.O) ,OR, ( ZT , LT, X. AND. H. 8T , 0 . 0 ) ) Na-H 
IF  (NT.LE,400)GOT0100$  PRINT  l.NTS  STOP 
F0RmAT(22H  run  error.  MERSON.  Na.IlO) 

XSaX  S  DO  no  Jal.NT  $  Q(J)»Y(J) 

CONTINUE 

HSaH*  OaX.H-ZT$  REa.TRUE, 

IF ( (Q.LT.O.O.AND.H.GE.O.O) .OR, (Q,aT,0,0.ANO,H,LE.0,0) )  GO  TO  210 
H=ZT-XS  flR=, FALSE. 

H3aH/3.0  $  DO  SIO  IS'«f*l,5  5  CALL  FUNC  <  NT .  X ,  Y .  F )  $  00  450  lal.NT 
QsH3*F(l)$  GOTO (30 1.30 2, 30 3. 304. 305) ,isw 
T ( I ) aRaQS  goto  400 
R=0,5* (O.T ( I > ) S  GOTO  400 

S(  1 )  =Ra3,0*(3S  Ra0,375*  (RoT  (I)  )  S  GO  TO  400 
T  {  I )  aRsT  (I  )  ♦4.0*0*  R«1  .5*  (R-SU  )  )  *  80  TO  400 
Ra0.5*(Q.T(l))S  O*AaS(2.0*R-1.5* (Q*S(I) ) ) 

Y(I)aG(I)*R  *  IFdSW.NE.S)  GO  TO  450  $  1 F  (.  NOT ,  BC )  GOTO  450  *  RaF5 
IF ( AflS  < Y ( I ) )  .GE. 0. 001 ) RaR*A8S( Y (I)  ) 

IF( (Q.LT.R) .0R..N0T.BX)G0T0  420  »  BRa.TRUE.S  BRa, FALSE,*  Ha0.5*H 
IF (ABS (H) .GE.HMI )GOTO  410  S  HaS I GNl ( H) *MMI i  flXa, FALSE. 

410  DO  411  Jal.NT  *  V(J)aG(J) 

411  CONTINUE  $  Kaxs*  GOTO  200 
4?0  IE  (Q. GE, 0. 031 ?5*R)BE», false. 

4S0  CONTINUE  *  GOTO t 50 1 . 5 1 0 . 50 3 ♦ 50 4 , 5 1 0 ) , I SW 
501  K«X+H3S  GOTO  510 

503  x«X*0.5*H3*  goto  510 

504  X«X.0.5*H 

510  CONTINUE  $  IEI.nOT.HC)  GO  TO  521 

IF ( .NOT, (BE. ANn.HH.ANO.BR) )  30  TO  ‘20  *  Hb2,0*M  »  flXa.TflUE, 

520  BHa.TRUE. 

521  IE (PR)  GO  TO  100  *  H«HS*  RETURN  *  END 
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SUBROUTINE  MATTNV 

OBTAINED  iPRON  COMPUTER  SUPPORT  DIVISION 
ABERDEEN  RESEARCH  AND  DEVELOPPENT  CENTER 


SUBROUTINE  NAT INV ( A, N«C«NNA8 »K,OET 1 

♦♦♦**  1 

OIPENSION  AINNAX, l),C( 1) 

MATINV  2 

NN  «  N 

MATINV  3 

KK  «  K 

MATINV  A 

IF  U-KKI  3il,l 

MATINV  5 

I 

N3  «  NN 

HAT INV  6 

IF  IKK)  2.4.2 

MATINV  7 

2 

ASSIGN  9  TO  NS 

MATINV  fa 

ASSIGN  13  TO  NT 

MATINV  9 

GOTO  5 

mat INVIO 

3 

N3  *  KK  «  NN  -  1 

MAT INVl  1 

A 

ASSIGN  1C  TO  NS 

MATINV12 

ASSIGN  lA  TO  NT 

MATINV  1  3 

5 

OET  »  l.O 

MAT  INV  I  *. 

00  IS  I  -  l.NN 

MAT  INVU 

IF  I  AM.  I  )  )  7.6.7 

HAT INV 16 

6 

hRITEIA. IT) 

MAT IMVl 7 

OET  -  0.0 

MATINV18 

GOTO  16 

HATINV19 

7 

T1  «  l.O/AI I. I ) 

MAT 1NV20 

CET  -  OET*AI I. I ) 

MATINV?! 

A(  1.  I  )  >  UO 

MAT INV22 

CO  8  J  >  1.N3 

HAT  I NV23 

A(  I.J)  -  A( I.J)«T1 

MAT  1 NV?4 

8 

CUNT  I NUE 

MAT  I NV25 

GOTO  N5.  (9,10) 

MAT  I NV2fa 

V 

C(  I )  -  Cl  I )*T1 

MAT  INV27 

10 

00  14  J  -  l.NN 

MAT INV?8 

IF  (I-J)  11,14.11 

MAT  INV?*.' 

1  1 

T1  ■  AU.  I  ) 

MAT INV 30 

A( J,  I  )  -  0.0 

MAT  I NV3  1 

CO  12  L  -  1,N3 

MAT  I NV3? 

AIJ.L)  -  AU.L)  -  T1*A(I,L) 

MAT  I  NV3  3 

12 

CONTINUE 

MAT  INV3<i. 

GOTO  NT,  I  13,14) 

MAT  I NV35 

1  3 

CIJ)  -  CIJ)  -  T1*CII) 

MAT  IN. '36 

l<i 

CONTINUE 

MAT INV37 

15 

CONT INUE 

MAT  I  N;'3R 

16 

RETURN 

HAT  I NV3  9 

17 

FC  NAT  1  16H  SINGULAR  MATRIX) 

MAT  I  NV'.C 

END 

MAT  I NVA  1 

'.)1 


No .  of 
Copies 


DISTRIBUTION  LIST 


Organization 
12  Conunander 

Defense  Technical  Info  Center 
ATTN:  DDC-DDA 
Cameron  Station 
Alexandria,  VA  22314 

1  Commander 

US  Army  Materiel  Development 
and  Readiness  Command 
AITN:  DRCDMD-ST 
5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

2  Commander 

US  Army  Armament  Research  and 
Development  Command 
ArrN;  DRDAR-TSS  (2  cys) 
Dover,  NJ  07801 

4  Commander 

US  Army  Armament  Research  and 
Development  Comirand 
ATTN:  DRDAR-LCA, 

Mr.  W.  R.  Benson 
DRDAR-LCA-F, 

Mr .  A .  Lo  eb 
DRDAR-I.CA-FA, 

Mr.  E.  Friedman 
Mr.  D.  Mertz 
Dover,  N.T  07801 

2  Commander 

US  Army  Armament  Research  .and 
Development  Command 
ATTN:  DRDAR-LCN, 

Mr.  F.  Saxe 
Mr.  F.  Seer bo 
Dover,  NJ  07801 

1  Comm-ander 

US  Aimy  Armament  Research  and 
Development  Command 
ATI'N:  imDAR-LCU,  Mr.  A.  Moss 
Dover,  N.i  07  801 


No.  of 

Copies  Organization 

1  Commander 

US  Army  Armament  Materiel 
Readiness  Command 
ATTN;  DRSAR-LEP-L,  Tech  Lib 
Rock  Island,  IL  61299 

1  Director 

US  Army  Armament  Research  and 
Development  Command 
Benet  Weapons  Laboratory 
ATTN:  DRDAR-LCB-TL 
Watervliet,  NY  12189 

1  Commander 

US  Army  Aviation  Research 
and  Development  Command 
ATTN:  DRSAV-R 
P.O.  Box  209 
St.  Louis,  MO  61366 

1  Director 

US  Army  Air  Mobility  Research 
and  Development  Laboratory 
Ames  Research  Center 
Moffett  Field,  CA  9403S 

1  Commander 

US  Army  Communications  Research 
and  Development  Command 
ATrN:  DRDCO-I>PA-SA 
Fort  Monmouth,  N.J  07703 

1  Coiamander 

US  Army  electronics  Research 
and  Devel  o])ment  Command 
Technical  Support  Activity 
ATFN:  DF.LSD-I. 

Fort  Monmoutii,  NJ  07703 

1  Commander 

US  Army  .Missile  Command 

ATFN:  DR  SMI -R 

Redstone  Ar.senal,  AI,  3.8800 
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DISTRIBUTION  LIST  (continued) 


No.  of  No.  of 

Copies  Organization  Copies  Organization 


1  Commander 

US  Army  Missile  Command 
ATTN:  DRSMI-YDL 
Redstone  Arsenal,  AL  35809 

1  Commander 

US  Army  Tank  Automotive 
Research  and  Development 
Command 

A'n'N:  DRDTA-UL 
Warren,  MI  48090 

1  Director 

US  Army  TRADOC  Systems 
Analysis  Activity 
ATTN:  ATAA-SL,  Tech  Lib 
V.'hite  Sands  Missile  Range 
NM  88002 

1  Commander 

US  Army  Yuma  Proving  Ground 
ATTN:  STF.YP-TMW, 

Mr.  W.  T.  Vomocil 
Yuma,  AZ  85.164 

1  Commander 

US  Arm)’  Research  Office 
ATTN:  CRD-AA-ldl 

1.0.  Box  12211 
Resi’.ireh  Tri.inj'Je  Park 
Nt:  27  70'.) 

1  Al-ATL  (Tech  l.ih) 

iigliii  Al-B,  i-l,  ,32542 

!  Afl'DL 

Wright-l’at  terson  APB, 
on  45433 

4  1)  i  r  ec  t  o  i' 

National  Aeronautics  and 
Space  Administration 
Ames  Research  enter 
ATI’N  :  Dr.  Car)’  Chapman 
Mr .  A .  So i  f  f 
Mr.  Murray  lobak 
lech  Lib 

Moffett  iield,  CA  9*1035 


2  Director 

Jet  Propulsion  Laboratory 
ATTN:  Tech  Lib 

Mr.  Peter  Joffe 
4800  Oak  Grove  Drive 
Pasadena,  CA  91103 

1  Commander 

David  W.  Taylor  Naval  Ship 
Research  and  Development 
Center 

ATTN;  Teen  Lib 
Bethesda,  MD  20084 

1  Commander 

Naval  Research  Laboratory 
ATTN:  Tech  Info  Div 
Washington,  D.C.  20375 

5  Commander 

Naval  Surface  Weapons  Center 
ATTN;  Dr.  Thomas  Clare 
Dr.  W.R.  Chadwick 
Dr.  W.G.  Soper 
Dr.  1'.  Moore 
Dr.  T.R.  Pejiitone 
Dahlgren,  VA  22448 

1  Commander 

Nax'al  Surface  Weapons  Center 
ATTN;  Code  7  30,  Tech  l.ih 
Silver  Spring,  MD  20910 

1  Commander 

Nav;il  Weapcais  Center 
A’Cl'N ;  Codi‘  233 
China  Lake,  CA  93555 

1  Arnold  Re  ns'ii  ch  Organ  i  c at  i  t)n  , 
1  nc  . 

I’rojeet  Support  and  Siu'cial 
Studies  Sect i on 
•Aoroch'nam )  c  s  Division 
Projects  Hraiich 
A ITN  ;  Dr,  .John  C.  Ad;iins  ,  .In. 
Arnold  M'S,  IN  3'’3S9 
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DISTRIBUTION  LIST  (continued) 


No.  of  No.  of 

Copies  Organization  Copies  Organization 


2  Director 

Sandia  Laboratories 
ATTN:  Division  13"! , 

Mr.  H.R.  Vaughn 
Mr.  A.E.  Hodapp,  Jr. 
Albuquerque,  NN'  87115 

1  Director 

National  Aercnautics  and 
Space  Admir LStration 
George  C.  Marshall  Space 
Flight  Center 
ATTN:  MS- I,  Library 
Huntsville,  AL  35812 

2  Director 

National  Aeronautics  and 
Space  Administration 
Langley  Research  Center 
ATTN:  MS- 185,  Tech  Lib 

Dr.  Clarence  Young 
Langley  Station 
Hampton,  VA  23365 

1  The  Analytic  Sciences 

Corporation  (TASC) 

ATTN;  Mr.  dames  1;.  Kai; 

6  .lacob  Way 
Reading,  M.A  01867 

1  i;.I.  elu  Font  de  .Nemours  and 

(Company 

Idiginoer  i  iig  Depai'tment 
A  I’TN  :  Di’.  Donald  W.  ^larciuaiait 
Wilmington  Dli  19898 

1  Exxon  Resea r'ch  Cientc-r 

AT'l'N :  Mr.  .lohii  St<  ven. 

Bo i 1 d i ng  28 
I’.d.  liox  ■';> 

Linden,  Nd  OtoSCj 

1  'denerai  I  leetrie  Comjnin)’ 

Armament  Systems  Department 
AT'l'N:  Ml'.  Robert  !!.  W^htte 

l.akesiile  Avenue 
Burl  i  ng!  on  ,  \  T  0.5  101 
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1  Oceanics,  Inc. 

ATTN:  Dr.  Theodore  R.  Goodman 
Technical  Industrial  Park 
Plainview,  NY  1IS03 

1  Schering  Corporation 

ATTN:  M,  Miller 
60  Orange  Street 
Bloomfield,  ,N.I  07003 

1  Systems  Technology,  Inc, 

ATTN:  Arlene  Muise,  Librarian 
13766  South  Hawthorne  Blvd. 
Hawthorne,  CA  90250 

1  Union  Research  Center 

ATTN:  Dr.  M.  A.  Selim 
P.O.  Box  76 
Brea,  CA  92621 

1  Case  Western  Reserve 

Uni versity 

ATT.N :  Professor  Philip  R 
Bevi ngton 

Cleveland,  OH  41106 

1  Oklahoma  State  llnivorsit>' 

Comjiuting  and  Information 
Sc  i  t'nces 

ATTN:  d,  P.  Cliandler 

Stillw.ater,  OK  74074 

1  Uni  vers  it  \’  of  Virg,  ini:i 

Di-partmeiit  •?!'  Idii;  i  ree i' i  n y. 

Sc i ence 

ATl'N:  I’rotessor  i  r;i  d,ico(>  -oii 

i'hornton  Hall 

(!har  '  ot  t  esv  i  1  1  e  ,  V A  .’2!)0  1 

-Abeialein  I'roving  (uiiiiml 
DT  r V  liS  VISA  \ 

A  I  IN.  |)k\.SV.|i 

DR.\  Si  -  ‘  li  ’  1 1  I  .  >j,i  1 1 

Cdi',  IISMliO'l 

\l IN  DRSil  iM  I 
I'l  I'.  USAC.'-I  ,  111  dr.  I  '..Mo,  I  \ 

M  IN:  ORO  M'  (II;  \ 
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USER  EVALUATION  OF  REPORT 


Please  take  a  few  minutes  to  answer  tlie  questions  liclow;  tear  out 
tills  sheet,  fold  as  indicated,  staple  or  tape  closed,  and  place 
in  tlie  mail.  Your  comments  will  provide  us  with  information  for 
imprevint;  future  reports. 

1 .  BRL  Report  Number  _ _ _  _ _ 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas,  etc.)  _  _  _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man-hours/contract  doll.ars  saved,  ojierating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate, 


.3.  Ceneral  Comments  (Indicate  what  ytm  tliink  should  In-  changed  to 
make  this  report  and  future  rejiort  s  of  this  type  irore  responsive 
to  your  need;;,  more  usable,  iiiijirove  readab  i  '  i  t  >• ,  etc., 


().  If  you  wouul  like  to  be  contacted  b\'  the  pi* r sonne  1  who  pri'pared 
this  report  to  raise  siiecifie  ipiestlvins  or  liiscuss  t  lu'  topic, 
please  fill  in  the  following  information. 

Name  : 


lelephone  Number: 


Hrgan i cat i on  Add  res  s ; 


